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We describe in detail how to implement a coarse-grained hybrid Molecular Dynamics and Stochas- 
tic Rotation Dynamics simulation technique that captures the combined effects of Brownian and 
hydrodynamic forces in colloidal suspensions. The importance of carefully tuning the simulation 
parameters to correctly resolve the multiple time and length-scales of this problem is emphasized. 
We systematically analyze how our coarse-graining scheme resolves dimensionless hydrodynamic 
numbers such as the Reynolds number Re, which indicates the importance of inertial effects, the 
Schmidt number Sc, which indicates whether momentum transport is liquid-like or gas-like, the Mach 
number Ma, which measures compressibility effects, the Knudsen number Kn, which describes the 
importance of non-continuum molecular effects and the Peclet number Pe, which describes the rel- 
ative effects of convective and diffusive transport. With these dimensionless numbers in the correct 
regime the many Brownian and hydrodynamic time-scales can be telescoped together to maximize 
computational efficiency while still correctly resolving the physically relevant physical processes. 
We also show how to control a number of numerical artifacts, such as finite size effects and solvent 
induced attractive depletion interactions. When all these considerations are properly taken into 
account, the measured colloidal velocity auto-correlation functions and related self diffusion and 
friction coefficients compare quantitatively with theoretical calculations. By contrast, these cal- 
culations demonstrate that, notwithstanding its seductive simplicity, the basic Langevin equation 
does a remarkably poor job of capturing the decay rate of the velocity auto-correlation function in 
the colloidal regime, strongly underestimating it at short times and strongly overestimating it at 
long times. Finally, we discuss in detail how to map the parameters of our method onto physical 
systems, and from this extract more general lessons - keeping in mind that there is no such thing 
as a free lunch - that may be relevant for other coarse-graining schemes such as Lattice Boltzmann 
or Dissipative Particle Dynamics. 

PACS numbers: 05.40.-a,82.70.Dd,47.11.-j,47.20.Bp 



I. INTRODUCTION 

Calculating the non- equilibrium properties of colloidal 
suspensions is a highly non-trivial exercise because these 
depend both on the short-time thermal Brownian mo- 
tion, and the long-time hydrodynamic behavior of the 
solvent Moreover, the hydrodynamic interactions 

are of a many-body character, and cannot usually be de- 
composed into a pairwise sum of inter-colloid forces. 

The fundamental difficulty of fully including the de- 
tailed solvent dynamics in computer simulations be- 
comes apparent when considering the enormous time and 
length-scale differences between mesoscopic colloidal and 
microscopic solvent particles. For example, a typical col- 
loid of diameter l/im will displace on the order of 10^° wa- 
ter molecules! Furthermore, a molecular dynamics (MD) 
scheme for the solvent would need to resolve time-scales 
on the order of 10~^^ s to describe the inter-molecular 
forces, while a colloid of diameter l/im in water diffuses 
over its own diameter in about Is. 

Clearly, simulating even an extremely crude molecular 
model for the solvent on the time-scales of interest is com- 
pletely out of the question: some form of coarse-graining 
is necessary. The object of this paper is to describe in 



detail one such scheme. But before we do so, we first 
briefly discuss a subset of the wide variety of different 
simulation techniques that have been devised to describe 
the dynamics of colloidal suspensions. 

At the simplest level, the effects of the solvent can be 
taken into account through Brownian dynamics (BD) 0], 
which assumes that collisions with the solvent molecules 
induce a random displacement of the colloidal particle 
positions, as well as a local friction proportional to the 
their velocity. Although, due to its simplicity, Brownian 
dynamics is understandably very popular, it completely 
neglects momentum transport through the solvent - as 
described by the Navier Stokes equations - which leads to 
long-ranged hydrodynamic interactions (HI) between the 
suspended particles. These HI may fall off as slowly as 
l/r, and can qualitatively affect the dynamical behavior 
of the suspension P, Q . 

Beginning with the pioneering work of Ermak and Mc- 
Cammon jj], who added a simple representation of the 
Oseen tensor Q to their implementation of BD, many 
authors have applied computational approaches that in- 
clude the HI by an approximate analytical form. The 
most successful of these methods is Stokesian Dynam- 
ics which can take into account higher order terms 
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in a multipole expansion of the HI interactions. Al- 
though some more sophisticated recent implementations 
of Stokesian Dynamics have achieved an 0{N log N) scal- 
ing with the number of colloids this method is still 
relatively slow, and becomes difficult to implement in 
complex boundary conditions. 

Another way to solve for the HI is by direct numer- 
ical simulation (DNS) methods, where the solid parti- 
cles are described by explicit boundary conditions for the 
Navier-Stokes equations 0. These methods are better 
adapted to non-Brownian particles, but can still be ap- 
plied to understand the effects of HI on colloidal dynam- 
ics. The fluid particle dynamics (FPD) method of Tanaka 
and Araki Q, and a related method by Yamamoto et 
al. are two important recent approaches to solving 
the Navier Stokes equations that go a step beyond stan- 
dard DNS, and simplify the problem of boundary con- 
ditions by using a smoothed step function at the colloid 
surfaces. In principle Brownian motion can be added to 
these methods |lo| . 

The difficulty of including complex boundary condi- 
tions in DNS approaches has stimulated the use of lattice- 
gas based techniques to resolve the Navier Stokes equa- 
tions. These methods exploit the fact that only a few con- 
ditions, such as (local) energy and momentum conserva- 
tion, need to be satisfied to allow the correct (thermo) hy- 
drodynamics to emerge in the continuum limit. Greatly 
simplified particle collision rules, easily amenable to ef- 
ficient computer simulation, are therefore possible, and 
complex boundary conditions, such as those that char- 
acterize moving colloids, are easier to treat than in DNS 
methods. The most popular of these techniques is Lat- 
tice Boltzmann (LB) where a linearized and pre-averaged 
Boltzmann equation is discretized and solved on a lat- 
tice Ladd has pioneered the application of this 
method to solid-fluid suspensions This is il- 
lustrated schematically in Fig. The solid particles 
are modeled as hollow spheres, propagated with New- 
tonian dynamics, and coupled to the LB solvent through 
a bounce-back collision rule on the surface. The fluctu- 
ations that lead to Brownian motion were modeled by 
adding a stochastic term to the stress tensor. Recently 
this has been criticized and an improved method to in- 
clude Brownian noise that does not suffer from some lat- 
tice induced artifacts was suggested ^ number of 
other groups have recently derived alternative ways of 
coupling a colloidal particle to a LB fluid [isl Ha Il7l | , in 
part to simulate the dynamics of charged systems. 

The discrete nature of lattice based methods can also 
bring disadvantages, particularly when treating fluids 
with more than one length-scale. Dissipative Particle Dy- 
namics (DPD) [isl IT^ is a popular off- lattice alternative 
that includes hydrodynamics and Brownian fluctuations. 
It can be viewed as an extension of standard Newtonian 
MD techniques, but with two important innovations: 

1) soft potentials that allow large time-steps and rapid 
equilibration; 

2) a Galilean invariant thermostat that locally conserves 




FIG. 1: Schematic picture depicting how a colloidal 
dispersion can be coarse-grained by replacing the solvent 
with a particle-based hydrodynamics solver such as Lattice 
Boltzmann (LB), Dissipative Particle Dynamics (DPD), the 
Lowe- Anderson thermostat, or Stochastic Rotation Dynamics 
(SRD). Each method introduces an effective coarse-graining 
length-scale that is chosen to be smaller than those of the 
mesoscopic colloids but much larger than the natural length- 
scales of a microscopic solvent. By obeying local momentum 
conservation, all four methods reproduce Navier-Stokes hy- 
drodynamics on larger length scales. For LB thermal fluctua- 
tions must be added in separately, but these emerge naturally 
for the other three methods. In this paper we focus on SRD, 
but many of the lessons learned should also apply to other 
particle based coarse-graining methods. 



momentum and therefore generates the correct Navier 
Stokes hydrodynamics in the continuum limit. 

In fact, these two methodological advances can be sep- 
arated. Soft potentials such as those postulated for in- 
novation 1) may indeed arise from careful equilibrium 
coarse-graining of complex fluids [2(ll l2ll . How- 

ever, their proper statistical mechanical interpretation, 
even for equilibrium properties, is quite subtle. For ex- 
ample, a potential that generates the correct pair struc- 
ture will not normally reproduce the correct virial pres- 
sure due to so called representability problems When 
used in dynamical simulations, the correct application of 
effective potentials is even more difficult to properly de- 
rive. For polymer dynamics, for instance, uncrossabil- 
ity constraints must be re-introduced to prevent coarse- 
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grained polymers from passing through one another |25j | . 
Depletion interactions in a multi component solution also 
depend on the relative rates of depletant diffusion coef- 
ficients and particle flow velocities "2^, '27*1 . At present, 
therefore, the statistical mechanical origins of innovation 
1) of DPD, the use of soft potentials out of equilibrium, 
are at best obscure. We are convinced that a correct 
microscopic derivation of the coarse-grained DPD repre- 
sentation of the dynamics, if this can indeed be done, 
will show that the interpretation of such soft potentials 
depends on dynamic as well as static (phase-space) av- 
erages. Viewing the DPD particles as static "clumps" 
of underlying fluid is almost certainly incorrect. It may, 
in fact, be more fruitful to abandon simple analogies to 
the potential energy of a Hamiltonian system, and in- 
stead view the interactions as a kind of coarse-grained 
self-energy p^ . 

Innovation 2), on the other hand, can be put on firmer 
statistical mechanical footing, see e.g. (29j, and can be 
usefully employed to study the dynamics of complex sys - 
tems with other types of inter-particle interactions |3Clj| . 
The main advantage of the DPD thermostat is that, by 
preserving momentum conservation, the hydrodynamic 
interactions that intrinsically arise from microcanonical 
MD are preserved for calculations in the canonical en- 
semble. Other thermostats typically screen the hydro- 
dynamic interactions beyond a certain length-scale [3lj |. 
For weak damping this may not be a problem, but for 
strong damping it could be. 

Simulating only the colloids with DPD ignores the 
dominant solvent hydrodynamics. While the solvent 
could be treated directly by DPD, as suggested in Fig.Q] 
(see also 113), this method is quite computationally ex- 
pensive because the "solvent" particles interact through 
pairwise potentials. 

In an important paper js^ , Lowe took these ideas fur- 
ther and combined the local momentum conservation of 
the DPD thermostat with the stochastic nature of the 
Anderson thermostat to derive a coarse-grained scheme, 
now called the Lowe- Anderson thermostat, which is much 
more efficient than treating the solvent with full DPD. It 
has recently been applied, for example, to polymer dy- 
namics |33 |. 

Independently, Malevanets and Kapral [s^ derived a 
method now called stochastic rotation dynamics (SRD) 
or also multiple particle collision dynamics (we choose 
the former nomenclature here). In many ways it resem- 
bles the Lowe- Anderson thermostat |33, or the much 
older direct simulation Monte Carlo (DSMC) method of 
Bird j3^,37||]. For all three of these methods, the parti- 
cles are ideal, and move in a continuous space, subject 
to Newton's laws of motion. At discrete time-steps, a 
coarse-grained collision step allows particles to exchange 
momentum. In SRD, space is partitioned into a rectan- 
gular grid, and at discrete time-steps the particles inside 
each cell exchange momentum by rotating their velocity 
vectors relative to the center of mass velocity of the cell 
(see Fig.nj. Similarly, in Lowe's method, a particle can. 



with a certain probability, exchange its relative velocity 
with another that lies within a certain radius. One can 
imagine other collision rules as well. As long as they lo- 
cally conserve momentum and energy, just as the lattice 
gas methods do, they generate the correct Navier Stokes 
hydrodynamics, although for SRD it is necessary to in- 
clude a grid-shift procedure to enforce Galilean invari- 
ance, something first pointed out by Ihle and KroU 3^. 
An important advantage of SRD is that its simplified 
dynamics have allowed the analytic calculation of several 
transport coefficients |3^ 113, |41| , greatly facifitating its 
use. In the rest of this paper we will concentrate on the 
SRD method, although many of our general conclusions 
and derivations should be easily extendible to the Lowe- 
Anderson thermostat and related methods summarized 
in Fig.Hl 

SRD can be applied to directly simulate fiow, as done 
for example in refs [i^ l43l |. but its stochastic nature 
means that a noise average must be performed to cal- 
culate flow lines, and this may make it less efficient 
than pre-averaged methods like LB. Where SRD becomes 
more attractive is for the simulation of complex particles 
embedded in a solvent. This is because in addition to 
long-ranged HI, it also naturally contains Brownian fiuc- 
tuations, and both typically need to be resolved for a 
proper statistical mechanical treatment of the dynamics 
of mesoscopic suspended particles. For example, SRD 
has been used to study polymers under flow |4J, [i^ |4g , 
the effect of hydrodynamics on protein folding and poly- 
mer collapse |47| , and the conformations of vesicles under 
flow |48l |. In each of the previously mentioned examples, 
the suspended particles are coupled to the solvent by par- 
ticipating in the collision step. 

A less coarse-grained coupling can be achieved by al- 
lowing direct collisions, obeying Newton's laws of mo- 
tion, between the SRD solvent and the suspended parti- 
cles. Such an approach is important for systems, such as 
colloidal suspensions, where the solvent and colloid time 
and length-scales need to be clearly separated. Male- 
vanets and Kapral derived such a hybrid algorithm 
that combined a full MD scheme of the solute-solute and 
solute-solvent interactions, while treating the solvent- 
solvent interactions via SRD. Early applications were to 
a two-dimensional many-particle swtem 50] , and to the 
aggregation of colloidal particles [51| . 

We have recently extended this approach, and applied 
it to the sedimentation of up to 800 hard sphere (HS) like 
colloids as a function of volume fraction, for a number of 
different values of the Peclet number js^- To achieve 
these simulations we adapted the method of Malevanets 
and Kapral |49| in a number of ways that were only briefly 
described in ref. [s^ . In the current paper we provide 
more in depth analysis of the simulation method we used, 
and describe some potential pitfalls. In particular we fo- 
cus on the different time and length-scales that arise from 
our coarse-graining approach, as well as the role of dimen- 
sionless hydrodynamic numbers that express the relative 
importance of competing physical phenomena. Very re- 
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cently, another similar study of colloidal sedimentation 
and aggregation has been carried out [S^l- Some of our 
results and analysis are similar, but some are different, 
and we will mention these where appropriate. 

After the introductory overview above, we begin in sec- 
tion ^] by carefully describing the properties of a pure 
SRD fluid, focusing on simple derivations that highlight 
the dominant physics involved for each transport coeffi- 
cient. Section lnil cxplains how to implement the coupling 
of a colloidal particle to an SRD solvent bath, and shows 
how to avoid spurious depletion interactions and how to 
understand lubrication forces. In Section llVI we analyze 
how optimizing the efficiency of particle based coarse- 
graining schemes affects different dimensionless hydrody- 
namic numbers, such as the Schmidt number, the Mach 
number, the Reynolds number, the Knudscn number, and 
the Peclet number. Section IVII describes the hierarchy 
of time-scales that determine the physics of a colloidal 
suspension, and compares these to the compressed hier- 
archy in our coarse-grained method. In section IVIII we 
tackle the question of how to map from a coarse-grained 
simulation, optimized for computational efficiency, to a 
real colloidal system. Finally after a Conclusions section, 
we include a few appendices that discuss amongst other 
things how to thermostat the SRD system (Appendix I); 
why the popular Langevin equation without memory ef- 
fects does a remarkably poor job of capturing both the 
long and the short time dynamics of the colloidal veloc- 
ity auto-correlation function (Appendix II); and, finally, 
how various physical properties and dimensionless num- 
bers scale for an SRD simulation, and how these compare 
to a colloid of radius 10 nm or 1 /^m in II2O (Appendix 
III). 



II. PROPERTIES OF A PURE SRD SOLVENT 

In SRD, the solvent is represented by a large number 
Nf of particles of mass mf. Here and in the following, 
we will call these "fluid" particles, with the caveat that, 
however tempting, they should not be viewed as some 
kind of composite particles or clusters made up of the 
underlying (molecular) fluid. Instead, SRD should be in- 
terpreted as a Navier Stokes solver that includes thermal 
noise. The particles are merely a convenient computa- 
tional device to facilitate the coarse-graining of the fluid 
properties. 



A. Simulation method for a pure solvent 

In the first (propagation) step of the algorithm, the 
positions and velocities of the fluid particles are propa- 
gated for a time Ate (the time between collision steps) 



by accurately integrating Newton's equations of motion. 



ruf 



dvj 
df 

dt 



(1) 
(2) 



Ti and Vi are the position and velocity of fluid particle i, 
respectively while is the total (external) force on par- 
ticle i, which may come from an external field such as 
gravity, or fixed boundary conditions such as hard walls, 
or moving boundary conditions such as suspended col- 
loids. The direct forces between pairs of fluid particles 
are, however, neglected in the propagation step. Herein 
lies the main advantage - the origin of the efficiency - 
of SRD. Instead of directly treating the interactions be- 
tween the fluid particles, a coarse-grained collision step 
is performed at each time-step At^. First, space is par- 
titioned into cubic cells of volume ag. Next, for each 
cell, the particle velocities relative to the center of mass 
velocity Vcm of the cell are rotated: 



R.(Vi - Vcm) ■ 



(3) 



R is a rotation matrix which rotates velocities by a fixed 
angle a around a randomly oriented axis. The aim of 
the collision step is to transfer momentum between the 
fluid particles while conserving the total momentum and 
energy of each cell. Both the collision and the streaming 
step conserve phase-space volume, and it has been shown 
that the single particle velocity distribution evolves to a 
Maxwell-Boltzmann distribution 35]. 

The rotation procedure can thus be viewed as a coarse- 
graining of particle collisions over space and time. Be- 
cause mass, momentum, and energy are conserved lo- 
cally, the correct hydrodynamic (Navier Stokes) equa- 
tions are captured in the continuum limit, including the 
effect of thermal noise [s^ . 

At low temperatures or small collision times Ate, 
the transport coefflcients of SRD, as originally formu- 
lated js^, show anomalies caused by the fact that fluid 
particles in a given cell can remain in that cell and partic- 
ipate in several collision steps |38j |. Under these circum- 
stances the assumption of molecular chaos and Galilean 
invariance are incorrect. However, this anomaly can be 
cured by applying a random shift of the cell coordinates 
before the collision step jsl, ■ 

It should also be noted that the collision step in SRD 
does not locally conserve angular momentum. As a con- 
sequence, the stress tensor a is not, in general, a symmet- 
ric function of the derivatives of the flow field (although 
it is still rotationally symmetric) [4l| . The asymmetric 
part can be interpreted as a viscous stress associated with 
the vorticity V x v of the velocity field v(r). The stress 
tensor will therefore depend on the amount of vorticity 
in the flow field. However, the total force on a fiuid el- 
ement is determined not by the stress tensor itself, but 
by its divergence V • a, which is what enters the Navier 
Stokes equations. Taking the divergence causes the ex- 
plicit vorticity dependence to drop out (the gradient of a 
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curl is zero). In principle, the stress tensor could be made 
symmetric by applying random rotations as well as ran- 
dom translations to the cell. Or, alternatively, angular 
momentum can be explicitly conserved by dynamically 
adapting the collision rule, as done by Ryder et al (54l |. 
who found no significant differences in fluid properties 
(although there is, of course, less flexibility in choosing 
simulation parameters). This result has been confirmed 
by some recent theoretical calculations [5^ , that demon- 
strate that the asymmetry of the stress tensor has only a 
few consequences, such as a correction to the sound wave 
attenuation associated with viscous dissipation of longi- 
tudinal density waves. These are not important for the 
fluid properties we are trying to model, and so, on bal- 
ance, we chose not to implement possible fixes to improve 
on angular momentum conservation. 



B. Transport coefficients and tlie dimensionless 
mean-free patli 

The simplicity of SRD collisions has facilitated the an- 
alytical calculation of many transport coefficients [s^ 
liol I41I l55l |. These analytical expressions are particu- 
larly useful because they enable us to efficiently tune the 
viscosity and other properties of the fluid, without the 
need for trial and error simulations. In this section we 
will summarize a number of these transport coefficients, 
where possible giving a simple derivation of the dominant 
physics. 



Units and the dimensionless mean-free path 



In this paper we will use the following units: lengths 
will be in units of cell-size oq, energies in units of ksT and 
masses in units of to/ (This corresponds to setting ao — I, 
ksT = 1 and mf = 1). Time, for example, is expressed 
in units of to = oq \fmjjk^ , the number density = 
7/ag and other derived units can be found in tableU We 
find it instructive to express the transport coefficients 
and other parameters of the SRD fluid in terms of the 
dimensionless mean-free path 



TABLE I: Units and simulation parameters for an SRD fluid. 
The parameters listed in the table all need to be fixed inde- 
pendently to determine a simulation. 



Basic Units 
ao = length 
ksT — energy 
mf — mass 



Derived Units 



to — ao 



Do = ^ = ao 



2 

ao 



ao 



ksT 
mf 

ksT _ 
mf 



time 



diflPusion constant 



kinematic viscosity 



7m f \/ mfkuT 

rjo = = 5 ~ viscosity 

toao ai 



Independent fluid simulation parameters 
7 — average number of particles per cell 
Ate ~ SRD collision time step 
a — SRD rotation angle 
L — box length 

Independent colloid simulation parameters 
AtuD ~ MD time step 

(Jcc ~ colloid-colloid collision diameter Eq. (Ilfcil 
ecc ~ colloid- colloid energy scale Eq. 1161 
(Jc/ — colloid-fluid collision diameter Eq. 1171 
ecf — colloid-fluid energy scale Eq. 1171 
Nc — number of colloids 
Mc — colloid mass 



2. Fluid self-diffusion constant 



a = ^Jm:.^, (4) 

ao y mf to 

which provides a measure of the average fraction of a cell 
size that a fluid particle travels between collisions. 

This particular choice of units helps highlight the basic 
physics of the coarse-graining method. The (nontrivial) 
question of how to map them on to the units of real 
physical system will be discussed in section IVIII 



A simple back of the envelope estimate of the self- 
diffusion constant £)y of a fluid particle can be obtained 
from a random- walk picture. In a unit of time to, a 
particle will experience 1/A collisions, in between which 
it moves an average distance Xao- Similarly, a heav- 
ier (tagged) particle of mass Mt, which exchanges mo- 
mentum with the fluid by participating in the coarse- 
grained collision step, will move an average distance 
Im — ^tcy/ksT/Mt = Xao/^/Mt between collisions. By 
viewing this motion as a random walk of step-size Im, 
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the diffusion coefficient follows: 

Do ^ Mt ' 



(5) 



expressed in units of Dq — ao/tQ — aoyJksT /mf. The 
diffusion coefficient D f for a pure fluid particle of mass 
TO/ is therefore given hy Dj/Dq w A. 

A more systematic derivation of the diffusion coeffi- 
cient of a fluid particle, but still within a random col- 
lision approximation, results in the following expres- 
sion I33.l5al: 



^^A 
Do 



2(1 — cos(q;)) \7 



(7-1) 


1" 




^ 2 



(6) 



The dependence on 7 is weak. If, for example, we 
take a = tt/2, the value used in this paper, then 
lim^^oo Df /Do = A, the same as Eq.(|Sll. 

A similar expression can be derived for the self- 
diffusion coefficient of a heavier tagged particle of mass 

Mt MM - 



£1 
Do 



Am/ 



2(1 - cos(a)) 



7 



(7) 



Note that this equation does not quite reduce to Eq. © 
for Mt/mf = 1 (because of slightly different approxi- 
mations), but the relative discrepancy decreases with in- 
creasing 7. 

While Eq. ^ is accurate for larger mean-free paths, 
where the random collision approximation is expected to 
be valid, it begins to show deviations from simulations for 
A ^ 0.6 IH^, when longer-time kinetic correlations begin 
to develop. RipoU et al JEB^ argue that these correlations 
induce interactions of a hydrodynamic nature that en- 
hance the diffusion coefficient for a fluid particle. For 
example, for A = 0.1 and a = jtt, they measured a fluid 
self-diffusion constant Df that is about 25% larger than 
the value found from Eq. ©. The enhancement is even 
more pronounced for heavier particles: for the same sim- 
ulation parameters they found that Dt was enhanced by 
about 75% over the prediction of Eq. ([7|) when Mt ^ 10. 

We note that coupling a large particle to the solvent 
through participation in the coarse-grained collision step 
leads to a diffusion coefficient which scales with mass as 
Dt/Do ~ i/Mt, whereas if one couples a colloid of ra- 
dius Rc to the solvent through direct MD collisions, one 
expects Dt/Do <io/Rc ~ {fnf/Mt)^. Moreover, it has 
been shown '50| that the effective hydrodynamic particle 
radius is approximately given by a w 2ao/(7r7). For any 
reasonable average number of fluid particles per cell, the 
effective hydrodynamic radius a of the heavier particle is 
therefore much less than the collision cell size ao ■ On the 
other hand, the hydrodynamic field is accurately resolved 
down to a scale comparable to ao {vide infra). What this 
implies is that this coupling method will yield correct 
hydrodynamic interactions only at large distances, when 
the colloids are more than several hydrodynamic radii 



apart. All the above suggests that some care must be 
taken when interpreting the dynamics of heavier parti- 
cles that couple through the coarse-grained collision step, 
especially when two or more heavy particles are in close 
proximity. 



3. Kinematic viscosity 

The spread of a velocity fluctuation dyjr) 
can be described by a diffusion equation |57| : 



dS\-{r) 
dt 



in a fluid 



(8) 



where v is the kinematic viscosity, which determines the 
rate at which momentum or vorticity "diffuses away" . 
The units of kinematic viscosity are vq = o,Q/to = 
aoy^kBT/ruf which are the same as those for particle 
self diffusion i.e. Do = vq- 

Momentum is transported through two mechanisms: 

1) By particles streaming between collision steps, lead- 
ing to a "kinetic" contribution to the kinematic viscos- 
ity J^kin- Since for this gas- like contribution the momen- 
tum is transported by particle motion, we expect Vkin 
to scale like the particle self-diffusion coefficient Df, i.e. 

Vkin/vo ^ A. 

2) By momentum being re-distributed among the parti- 
cles of each cell during the collision step, resulting in a 
"collisional" contribution to the kinematic viscosity i^coi- 
This mimics the way momentum is transferred due to 
intcr-particle collisions, and would be the dominant con- 
tribution in a dense fluid such as water at standard tem- 
perature and pressure. Again a simple random-walk ar- 
gument explains the expected scaling in SRD: Each col- 
lision step distributes momentum among particles in a 
cell, making a step-size that scales like cq. Since there are 
1/A collision steps per unit time to, this suggests that the 
collisional contribution to the kinematic viscosity should 
scale as Vcoi/vq ~ 1/A. 

Accurate analytical expressions for the kinematic vis- 
cosity V = Vkin + Vcoi of SRD have been derived "4^, '55| , 
and these can be rewritten in the following dimensionless 
form: 



=0 /kin (7, a) 

vo 3 

— = /col (7, a). 



(9) 
(10) 



where the dependence on the collisional angle a and fluid 
number density 7 is subsumed in the following two fac- 
tors: 



/kin (7,") 



157 



(7- 1 + e-'>')(4- 2cos(a) - 2cos(2a)) 2' 

(11) 



a) = (1 - cos(a))(l - 1/7 + e-V7)- 



(12) 



7 



These factors only depend weakly on 7 for the typical 
parameters used in simulations. For example, at a = 
7r/2 and 7 = 5, the angle and number density we use in 
this paper, they take the values /]J'j,j(5, 7r/2) = 1.620 and 
fco\{5, 7r/2) — 0.801. For this choice of collision angle a, 
they monotonically converge to /kin = /coi = 1 in the 
limit of large 7. 



and the density 7 can be varied independently. More- 
over, the coUisional contribution to the viscosity adds a 
new dimension, allowing a much wider range of physical 
fluids to be modeled than just simple gases. 



III. COLLOID SIMULATION METHOD 



4- Shear viscosity 

Suppose the fluid is sheared in the x direction, with the 
gradient of the average flow field in the y direction. Two 
neighboring fluid elements with different y-coordinates 
will then experience a friction force in the x-direction, as 
expressed by the xy component of the stress tensor a, 
which for a simple liquid is linearly proportional to the 
instantaneous flow field gradient. 



rj- 



dy 



(13) 



The coefficient of proportionality r/ is called the shear 
viscosity, and is it related to the kinematic viscosity by 
77 — PfV, where pf = nif^/aQ is the fluid mass density. 
From Eas. IHHT^it follows that the two contributions to 
the shear viscosity can be written in dimensionless form 
as: 



Vkin 

Vo 

Vcol 

Vo 



A7 



/ki„(7,a) 



Y^/c';i(7,«). 



(14) 
(15) 



where rjo = m/aoto — y^m^^fc^jT/ag is the unit of shear 
viscosity. 

In contrast to the expressions for the diffusion of a 
fluid particle or a tagged particle, Eqs. © - l|15(l com- 
pare quantitatively to simulations over a wide range of 
parameters [39l l40l l5a| . For the parameters we used in 
our simulations in |52L IH^ . i.e. A = 0.1, a = tt/2,j = 5, 
the collisional contribution to the viscosity dominates; 
i^kin = 0. 0541^0 and z^coi = 0.45fo. This is typical for 
A <C 1, where ly and rj can be taken to a good first ap- 
proximation by the collisional contribution only. In fact 
throughout this paper we will mainly focus on this small 
A limit. 

It is also instructive to compare the expressions de- 
rived in this section to what one would expect for simple 
gases, where, as famously first derived and demonstrated 
experimentally by Maxwell in the 1860's 59], the shear 
viscosity is independent of density. This result differs 
from the kinetic (gas-like) contribution to the viscosity 
in Eq. H14I) . because in a real dilute gas the mean- free 
path scales as A oc I/7, canceling the dominant density 
dependence in rjkin oc 7A. The same argument explains 
why the self-diffusion and kinematic viscosity of the SRD 
fluid are, to first order, independent of 7, while in a gas 
they would scale as I/7. In SRD, the mean-free path A 



Malevanets and Kapral l4!l| flrst showed how to imple- 
ment a hybrid MD scheme that couples a set of colloids 
to a bath of SRD particles. In this section, we expand on 
their method, describing in detail the implementation we 
used in ref. 52J. We restrict ourselves to HS like colloids 
with steep interparticle repulsions, although attractions 
between colloids can easily be added on. The colloid- 
colloid and colloid-fluid interactions, (fcc{f) and ipcf{r) 
respectively, are integrated via a normal MD procedure, 
while the fluid-fluid interactions are coarse-grained with 
SRD. Because the number of fluid particles vastly out- 
numbers the number of HS colloids, treating their in- 
teractions approximately via SRD greatly speeds up the 
simulation. 



A. Colloid-colloid and colloid-solvent interactions 

Although it is possible to implement an event-driven 
dynamics of HS colloids in an SRD solvent '6^ , here we 
approximate pure HS colloids by steep repulsive interac- 
tions of the WCA form j^: 



(r < 2i/24ae 



^,,(r) = \ ^'^^ [[—) ~ [ — ) + 4 

I (r > 21/24^^^; 

(16) 

Similarly, the colloid-fluid interaction takes the WCA 
form: 



4ec/ 




(r < 2^/^a,j) 



(r > 2VV,;). 
(17) 

The mass my of a fluid particle is typically much smaller 
than the mass Ale of a colloid, so that the average ther- 
mal velocity of the fluid particles is larger than that 
of the colloid particles by a factor Mcjmf. For this 
reason the time-step A^md is usually restricted by the 
fluid-colloid interaction (|17|l . allowing fairly large ex- 
ponents n for the colloid-colloid interaction (/?cc('") — 
^^cc [{(Jcc/rT' - {acc/rY' + 1/4] . We choose n = 24, 
which makes the colloid-colloid potential steep, more like 
hard-spheres, while still soft enough to allow the time- 
step to be set by the colloid-solvent interaction. 

The positions and velocities of the colloidal spheres are 
propagated through the Velocity Verlet algorithm |62 
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with a time step AtMD^ 

R,(^ + AtMD) = (<) + V, (t) AiMD 
V,(t + A<MD) = V,(t) 



2M, 

F, (t) +F,{t + AtMo) 



2M, 



AijviD- 
(19) 



Ri and are the position and velocity of colloid i, re- 
spectively. Fi is the total force on that colloid, exerted by 
the fluid particles, an external field, such as gravity, ex- 
ternal potentials such as repulsive walls, as well as other 
colloids within the range of the interaction potential (|16|l . 

The positions and velocities of SRD particles are 
updated by similarly solving Newton's Eqns. (|1I2|I every 
time-step ISImd , and with the SRD procedure of Eq. Q 
every time-step Ate- 

Choosing both AtMD and Ate as large as possible en- 
hances the efficiency of a simulation. To first order, each 
timestep is determined by different physics AtM d by the 
steepness of the potentials, and Ate by the desired fluid 
properties, and so there is some freedom in choosing their 
relative values. We used Atc/AtMB = 4 for our simula- 
tions of sedimentation [52i] , but other authors have used 
ratios of 50 |49l or even an order of magnitude larger 
than that |53|. Later in the paper we will revisit this 
question, linking the time-steps to various dimensionless 
hydrodynamic numbers and Brownian time-scales. 



B. Stick and slip boundary conditions 

Because the surface of a colloid is never perfectly 
smooth, collisions with fluid particles transfer angular 
as well as linear momentum. As demonstrated in Fig. |21 
the exact molecular details of the colloid-fluid interac- 
tions may be very complex, and mediated via co- and 
counter-ions, grafted polymer brushes etc. . . . However, 
on the time and length-scales over which our hybrid MD- 
SRD method coarse-grains the solvent, these interactions 
can be approximated by stick boundary conditions: the 
tangential velocity of the fluid, relative to the surface of 
the colloid, is zero at the surface of the colloid [63L l63|. 
For most situations, this boundary condition should be 
sufficient, although in some cases, such as a non-wetting 
surface, large slip-lengths may occur [65| . 

In computer simulations, stick boundary conditions 
may be implemented by bounce-back rules, where both 
parallel and perpendicular components of the relative ve- 
locity are reversed upon a collision with a surface. These 
have been applied by Lamura et al [4^ . who needed to 
modify the bounce-back rules slightly to properly repro- 
duce stick boundaries for a Poiseuille flow geometry. 

Stick boundaries can also be modeled by a stochastic 
rule. After a collision, the relative tangential velocity 
Vt and relative normal velocity Vn are taken from the 




FIG. 2: Schematic picture depicting how a fluid molecule in- 
teracts with a coUoid, imparting both linear and angular mo- 
mentum. Near the colloidal surface, here represented by the 
shaded region, there may be a steric stabilization layer, or a 
double-layer made up of co- and counter-ions. In SRD, the 
detailed manner in which a fluid particle interacts with this 
boundary layer is represented by a coarse-grained stick or slip 
boundary condition. 



distributions: 



P{Vn) OC Vnexp{-f3vl) 

P{vt) OC exp {-I3v1) , 



(20) 
(21) 



so that the colloid acts as an additional thermostat [63 ■ 
Such stochastic boundary conditions have been used for 
colloidal particles by Inoue et al. [H^ and Hecht et al. [s^ ■ 
We have systematically studied several implementations 
of stick boundary conditions for spherical colloids ^5^ . 
and derived a version of the stochastic boundary condi- 
tions which reproduces linear and angular velocity corre- 
lation functions that agree with Enskog theory for short 
times, and hydrodynamic mode-coupling theory for long 
times. We argue that the stochastic rule of Eq. H20|l is 
more like a real physical colloid - where fluid-surface 
interactions are mediated by steric stabilizing layers or 
local CO- and counter-ion concentrations - than bounce- 
back rules are [s^ . 

Nevertheless, in this paper, many examples will be for 
radial interactions such as those described in Eq. H17|l. 
These do not transfer angular momentum to a spherical 
colloid, and so induce effective slip boundary conditions. 
For many of the hydrodynamic effects we will discuss here 
the difference with stick boundary conditions is quanti- 
tative, not qualitative, and also well understood. 



C. Depletion and lubrication forces 

1. Spurious depletion forces induced by the fluid 

We would like to issue a warning that the additional 
fluid degrees of freedom may inadvertently introduce de- 
pletion forces between the colloids. Because the number 
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density of SRD particles is much higher than that of the 
colloids, even a small overlap between two colloids can 
lead to enormous attractions. 

For low colloid densities the equilibrium depletion in- 
teraction between any two colloids caused by the presence 
of the ideal fluid particles is given by [s^, Ell : 

$dcpl(c?) = njkBT [yexci(rf) - Kxci(oo)] , (22) 

where rif = 7/ag is the number density of fluid particles 
and Vexci{d) is the (free) volume excluded to the fluid by 
the presence of two colloids separated by a distance d. 
The latter is given by 

T4xci(d) = J d^r {1 - exp [~(3ipcf (r - ri) - (3ipcf (r - ra)]} 

(23) 

where |ri — ra] = d. An example is given in Fig. |21 where 
we have plotted the resulting depletion potential for the 
colloid-solvent interaction H17() . with Ccc = Cc/ = i.bksT 
as routinely used in our simulations, as well as the deple- 
tion interaction resulting from a truly HS colloid-solvent 
interaction. The latter can easily be calculated analyti- 
cally, with the result 

ford<2crc/. (24) 

For the pure HS interactions, one could take CTcc > 2(Tc/ 
and the depletion forces would have no effect. But for 
interactions such as those used in Eqs. l(T^ and ((T7jl . 
the softer repulsions mean that inter-colloid distances 
less than CTcc are regularly sampled. A more stringent 
criterion of CTcc must therefore be used to avoid spuri- 
ous depletion effects. As an example, we re-analyze the 
simulations of ref. 51], where a WCA form with n = 6 
was used for the fluid-colloid interactions, and a normal 
Lennard Jones n — 6 potential with (Tc/ > ^o^cc was used 
for the colloid-colloid interactions. The authors found 
differences between "vacuum" calculations without SRD 
particles, and a hybrid scheme coupling the colloids to an 
SRD solvent. These were correctly attributed to "solvent 
induced pressure" , which we quantify here as depletion 
interactions. For one set of their parameters, acf = ^(Jcc, 
the effect is mainly to soften the repulsion, but for their 
other parameter set: acf — 0.65f7cc, depletion attractions 
induce an effective attractive well-depth of over 40fcBT! 

Since the depletion potentials can be calculated an- 
alytically, one might try counteracting them by intro- 
ducing a compensating repulsive potential of the form: 
$comp = — $dcpi between the colloids. However, there 
are three problems with this approach: Firstly, at higher 
colloid packing fractions, three and higher order inter- 
actions may also need to be added, and these are very 
difflcult to calculate. Secondly, the depletion interactions 
are not instantaneous, and in fact only converge to their 
equilibrium average algebraically in time "5^ . While this 
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FIG. 3: (Color online) Effective depletion potentials induced 
between two colloids by the SRD fiuid particles, for HS fluid- 
colloid (solid line) and WCA fluid-colloid (dashed line) in- 
teractions taken from Eq. 1171 . The interparticle distance is 
measured in units of a^f- Whether these attractive poten- 
tials have important effects or not depends on the choice of 
diameter CTcc in the bare colloid- colloid potential 1161 . 

is not a problem for equilibrium properties, it will in- 
troduce errors for non-equilibrium properties. Finally, 
when external flelds drive the colloid, small but persis- 
tent anisotropics in the solvent density around a colloid 
may occur |23|. Although these density variations are 
(and should be) small, the resulting variations in deple- 
tion interactions can be large. 

To avoid these problems, we routinely choose the 
colloid-fluid interaction range acf slightly below half the 
colloid diameter acc/'^- More precisely, we ensure that 
the colloid-colloid interaction equals 2.5kBT at a distance 
d where the depletion interactions have become zero, i.e., 
at a distance of twice the colloid-solvent interaction cut- 
off radius. Smaller distances will consequently be rare, 
and adequately dealt with by the compensation poten- 
tial. This solution may be a more realistic representa- 
tion anyhow, since in practice for charge and even for 
sterically stabilized colloids, the effective colloid-colloid 
diameter acc is expected to be larger than twice the ef- 
fective colloid-fluid diameter acf - This is particularly so 
for charged colloids at large Debye screening lengths. 

2. Lubrication forces depend on surface details 

When two surfaces approach one another, they must 
displace the fluid between them, while if they move apart, 
fluid must flow into the space freed up between the sur- 
faces. At very short inter-surface distances, this results 
in so-called lubrication forces, which are repulsive for col- 
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loids approachingeach other, and attractive for colloids 
moving apart P, 0, E^l • These forces are expected to be 
particularly inipqrtant for driven dense colloidal suspen- 
sions; see e.g. j69| for a recent review. 

An additional advantage of our choice of diameters ad 
above is that more fluid particles will fit in the space 
between two colloids, and consequently the lubrication 
forces will be more accurately represented. It should be 
kept in mind that for colloids, the exact nature of the 
short-range lubrication forces will depend on physical de- 
tails of the surface, such as its roughness, or presence of 
a grafted polymeric stabihzing layer (70|, |7j, 0|. For 
perfectly smooth colloids, analytic limiting expressions 
for the lubrication forces can be derived [63j, showing a 
divergence at short distances. We have confirmed that 
SRD resolves these lubrication forces down to surpris- 
ingly low interparticle distances. But, at some point, this 
will break down, depending of course on the choice of sim- 
ulation parameters (such as tTc//ao, A, and 7), as well as 
the details of the particular typ e of colloidal particles that 
one wishes to model [23,0, IzS- An explicit analytic cor- 
rection could be applied to properly resolve these forces 
for very small distances, as was recently implemented for 
Lattice Boltzmann |73l] . However, in this paper, we will 
assume that our choice of ctc/ is small enough for SRD to 
sufficiently resolve lubrication forces. The lack of a com- 
plete divergence at very short distances may be a bet- 
ter model of what happens for real colloids anyway. For 
dense suspensions under strong shear, explicit lubrication 
force corrections as well as other short-ranged colloid- 
colloid interactions arising from surface details such as 
polymer coats will almost certainly need to be put in by 
hand, see e.g. ref. for further discussion of this subtle 
problem. 



D. Test of static properties 

The equilibrium properties of a statistical mechanical 
system should be independent of the detailed dynamics 
by which it explores phase space. Therefore one requisite 
condition imposed on our coarse-grained dynamics is that 
it reproduces the correct static properties of an ensemble 
of colloids. 

An obvious property to measure is the radial distribu- 
tion function g{r) In Fig. 0] we depict some g(r)s 
obtained from SRD simulations at a number of differ- 
ent densities, and compare these to similar simulations 
with standard Brownian dynamics. The colloids interact 
via potentials of the form \i<6\ - H17|) with cjcf ~ 2ao, 
(Tec — 4.3ao, and ecc = ^cf = 2.5kBT. For the SRD we 
used 7 = 5, a = ^TT, A^MD — 0-025to and Ate — O.l^o, 
implying A — 0.1. For the BD we used a background 
friction calculated from the effective hydrodynamic ra- 
dius {vide infra) and a time-step equal to the A^md used 
above. The colloids were placed in a box of dimensions 
32 X 32 X 96aQ, and the number of particles was varied 
from 64 to 540 in order to achieve different packing frac- 
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FIG. 4: (Color online) Colloid radial distribution functions 
g{r) for various colloid volume fractions (jjc. The simulations 
were performed with BD (black lines) and SRD (red lines) 
simulations, and are virtually indistinguishable within statis- 
tical errors. 



tions 

(l)c = ^TTUcale, (25) 

where the subscript in (j)c refers to volume fraction based 
on the colloid-colloid interaction, to distinguished it from 
the volume fraction (f) based on the colloid hydrodynamic 
radius fF3 |. From Fig. H is clear that the two simula- 
tions are indistinguishable within statistical errors, as 
required. Even though (Tc/ < ^CTcc, we still found it 
necessary to include an explicit compensating potential 
without which the g{r) generated by SRD is increased 
noticeably at contact when compared to BD simulations. 

Finally, we note that any simulation will show finite 
size effects. These may be thermodynamic as well as 
dynamic. If there are thermodynamic finite-size effects 
(for example due to a long correlation length caused by 
proximity to a critical point), then a correct simulation 
technique will show the same behaviour regardless of the 
underlying dynamics. 



IV. DIMENSIONLESS NUMBERS 

Different regimes of hydrodynamic behavior can be 
characterized by a series of dimensionless numbers that 
indicate the relative strengths of competing physical pro- 
cesses js^l- If two distinct physical systems can be de- 
scribed by the same set of hydrodynamic numbers, then 
their flow behavior will be governed by the same physics, 
even if their time and length scales differ by orders of 
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TABLE II: Dimensionless hydrodynamic numbers for col- 
loidal suspensions. 

„ , . , coUisional momentum transport „ v ^ — _ 
Schmidt — Sc = —— Eq. 11261 

Df 



than the fluid diffusion coefficient, i.e. Dcoi ^ -D/ to 
achieve the correct time-scale separation, as we will see 
in Section IvH 



Mach 



Reynolds 



Knudsen 



Peclet 



kinetic momentum transport 

flow velocity 
sound velocity 

inertial forces 



Ma 



viscous forces 

mean free path 
particle size 

convective transport 
diffusive transport 



Re = ^ Eq. m 



Kn 



Pe 



Eq. 



magnitude. In other words, a macroscopic boulder sedi- 
menting in viscous molten magma and a mesoscopic col- 
loid sedimenting in water may show similar behavior if 
they share key hydrodynamic dimensionless numbers. It 
is therefore very instructive to analyze where a system de- 
scribed by SRD sits in "hydrodynamic parameter space" . 
In this section we review this parameter space, one "di- 
mension" at a time, by exploring the hydrodynamic num- 
bers listed in Table ITTI 



Schmidt number 



The Schmidt number 



Sc 



V 

Id 



(26) 



is important for characterizing a pure fluid. It expresses 
the rate of diffusive momentum transfer, measured by 
the kinematic viscosity v, relative to the rate of diffusive 
mass transfer, measured by the fluid particle self-diffusion 
coefficient Df. For a gas, momentum transport is domi- 
nated by mass diffusion, so that Sc « 1, whereas for a liq- 
uid, momentum transport is dominated by inter-particle 
collisions, and Sc ^ 1. For low Schmidt numbers, the 
dynamics of SRD is indeed "gas-like" , while for larger 
Sc, the dynamics shows collective behavior reminiscent of 
the hydrodynamics of a liquid [s^ . This distinction will 
be particularly important for simulating the dynamics of 
embedded particles that couple to the solvent through 
participation in the collision step. For particles that 
couple directly, either through a potential, or through 
bounce-back or stochastic boundary conditions, this dis- 
tinction is less important because the fluid-particle dif- 
fusion coefficient Df does not directly enter the Navier 
Stokes equations. Of course other physical properties can 
be affected. For example, we expect that the self diffu- 
sion coefficient of a colloid Dr 



Eq. m 



1. Dependence of Schmidt number on simulation 
parameters 

From Eqs. © - H12() it follows that the Schmidt num- 
ber can be rewritten as 



Eq. (031 



Sc: 



1 



18A2' 



(27) 



where we have ignored the dependence on 7 and a since 
the dominant scaling is with the dimensionless mean-free 
path A. For larger A, where the kinetic contributions 
to the viscosity dominate, the Sc number is small, and 
the dynamics is gas-like. Because both the fluid diffu- 
sion coefficient © and the kinetic contribution to the 
viscosity (j^ scale linearly with A, the only way to ob- 
tain a large Sc number is in the limit A <C 1 where the 
coUisional contribution to ly dominates. 

Low Sc numbers may be a more general characteris- 
tic of particle-based coarse-graining methods. For com- 
putational reasons, the number of particles is normally 
greatly reduced compared to the solvent one is modeling. 
Thus the average inter-particle separation and mean-free 
path are substantially increased, typically leading to a 
smaller kinematic viscosity and a larger fluid diffusion 
coefficient Df. With DPD, for example, one typically 
finds Sc « 1 [2^ |33|| . It is therefore difficult to achieve 
large Sc numbers, as in a real liquid, without sacrific- 
ing computational efficiency. But this may not always 
be necessary. As long as momentum transport is clearly 
faster than mass transport, the solvent should behave in 
a liquid-like fashion. With this argument in mind, and 
also because the Sc number does not directly enter the 
Navier Stokes equations we want to solve, we use A = 0.1 
in this paper and in ref. [s^ . which leads to a Sc ~ 5 
for 7 = 5 and a = 7t/2. This should be sufficient for 
the problems we study. For other systems, such as those 
where the suspended particles are coupled to the solvent 
through the collision step, more care may be needed to 
ensure that the Sc number is indeed large enough (56j |. 



B. Mach number 

The Mach number measures the ratio 



Ma: 



(28) 



between Vg, the speed of solvent or colloid flow, and 
Cf — \J {b/?>){kBT /mf), the speed of sound. In contrast 
to the Schmidt number, which is an intrinsic property of 
the solvent, it depends directly on flow velocity. The 
Ma number measures compressibility effects j57| since 



12 



sound speed is related to the compressibility of a liq- 
uid. Because c/ in many liquids is of order 10"^ m/s, the 
Ma numbers for physical colloidal systems are extremely 
small under normally achievable flow conditions. Just 
as for the Sc number, however, particle based coarse- 
graining schemes drastically lower the Ma number. The 
particle mass is typically much greater than the mass 
of a molecule of the underlying fluid, resulting in lower 
velocities, and moreover, due to the lower density, colli- 
sions also occur less frequently. These effects mean that 
the speed of sound is much lower in a coarse-grained sys- 
tem than it is in the underlying physical fluid. Or, in 
other words, particle based coarse-graining systems are 
typically much more compressible than the solvents they 
model. 

Ma number effects typically scale with Ma^ j57,] , and so 
the Ma number does not need to be nearly as small as for 
a realistic fluid to still be in the correct regime of hydro- 
dynamic parameter space. This is convenient, because to 
lower the Ma number, one would need to integrate over 
longer fluid particle trajectories to allow, for example, 
a colloidal particle to flow over a given distance, making 
the simulation computationally more expensive. So there 
is a compromise between small Ma numbers and compu- 
tational efficiency. We limit our Ma numbers to values 
such that Ma < 0.1, but it might be possible, in some 
situations, to double or triple that limit without causing 
undue error. For example, incompressible hydrodynam- 
ics is used for aerodynamic flows up to such Ma numbers 
since the errors are expected to scale as 1/(1 — Ma^) 
When working in units of m/ and fc^T, the only way to 
keep the Ma number below our upper limit is to restrict 
the maximum flow velocity to Vs O.lc/. The flow veloc- 
ity itself is, of course, determined by the external flelds, 
but also by other parameters of the system. 



C. Reynolds number 

The Reynolds number is one of the most important di- 
mensionless numbers characterizing hydrodynamic flows. 
Mathematically, it measures the relative importance of 
the non- linear terms in the Navier-Stokes equations [s^ . 
Physically, it determines the relative importance of iner- 
tial over viscous forces, and can be expressed as: 



the same as the kinematic time 



-^ = Rets 

V 



Re 



(29) 



where a is a length-scale, in our case, the hydrodynamic 
radius of a colloid, i.e. a ~ CTc/, and Vs is a flow velocity. 

For a spherical particle in a flow, the following heuristic 
argument helps clarify the physics behind the Reynolds 
number: If the Stokes time 

ts = ^, (30) 

Vs 

it takes a particle to advect over its own radius is about 



(31) 



it takes momentum to diffuse over that distance, i.e. 
Ke—T,y/ts « 1, then the particle will feel vorticity ef- 
fects from its own motion a distance Uc/ away, leading to 
non-linear inertial contributions to its motion. Since hy- 
drodynamic interactions can decay as slowly as 1/r, their 
influence can be non-negligible. If, on the other hand. Re 
<C 1, then vorticity will have diffused away and the par- 
ticle will only feel very weak hydrodynamic effects from 
its own motion. 

Exactly when inertial flnite Re effects become signif- 
icant depends on the physical system under investiga- 
tion. For example, in a pipe, where the length-scale a in 
Eq. (|29|l is its diameter, the transition from simpler lami- 
nar to more complex turbulent flow is at a pipe Reynolds 
number of Re « 2000 [il,!?^. On the other hand, for a 
single spherical particle, a non-linear dependence of the 
friction ^ on the velocity Vg, induced by inertial effects, 
starts to become noticeable for a particle Reynolds num- 
ber of Re « 1 j57j . while deviations in the symmetry 
of the streamlines around a rotating sphere have been 
observed in calculations for Re « 0.1 |7q . 

For typical colloidal suspensions, where the particle 
diameter is on the order of a few /im down to a few 
nm, the particle Reynolds number is rarely more than 
10^'^. In this so-called Stokes regime viscous forces dom- 
inate and inertial effects can be completely ignored. The 
Navier-Stokes equations can be replaced by the linear 
Stokes equations so that analytic solutions are easier 
to obtain 63\. However, some of the resulting behavior is 
non-intuitive for those used to hydrodynamic effects on 
a macroscopic scale. For example, as famously explained 
by Purcell in a talk entitled "Life at low Reynolds num- 
bers" j77j, many simple processes in biology occur on 
small length scales, well into the Stokes regime. The 
conditions that bacteria, typically a few /im long, expe- 
rience in water are more akin to those humans would ex- 
perience in extremely thick molasses. Similarly, colloids, 
polymers, vesicles and other small suspended objects are 
all subject to the same physics, their motion dominated 
by viscous forces. For example, if for a colloid sediment- 
ing at 1 /Ltm/s, gravity were instantaneously turned off, 
then the Stokes equations suggest that it would come to a 
complete halt in a distance significantly less than one A, 
reflecting the irrelevance of inertial forces in this low Re 
number regime. It should be kept in mind that when the 
Stokes regime is reached because of small length scales (as 
opposed to very large viscosities such as those found for 
volcanic lava flows), then thermal fluctuations are also 
important. These will drive diffusive behavior [7^. In 
many ways SRD is ideally suited for this regime of small 
particles because the thermal fluctuations are naturally 
included [H. 
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1. Dependence of the Reynolds number on simulation 
parameters 

From Eqs. and it follows that the Reynolds 
number for a colloid of hydrodynamic radius a « acf can 
be written as: 



Re 



3 \ af) J \ v 



(32) 



Equating the hydrodynamic radius a to CTc/ from the 
fluid-colloid WCA interaction of Eq. (flTjl is not quite 
correct, as we will see in section but is a good enough 
approximation in light of the fact that these hydrody- 
namic numbers are qualitative rather than quantitative 
indicators. 

In order to keep the Reynolds number low, one must 
either use small particles, or very viscous fluids, or low 
velocities Vg. The latter condition is commensurate with 
a low Ma number, which is also desirable. For small 
enough A (and ignoring for simplicty the factor /^oi(7, a) 
in Eq. (llUfl ) the Re number then scales as 



Re: 



23 Ma ^ A. 

ao 



(33) 



Again, we see that a smaller mean- free path A, which 
enhances the coUisional viscosity, also helps bring the Re 
number down. This parameter choice is also consistent 
with a larger Sc number. Thus larger Sc and smaller Ma 
numbers both help keep the Re number low for SRD. 
In principle a large viscosity can also be obtained for 
large A, which enhances the kinetic viscosity, but this 
choice also lowers the Sc number and raises the Knudsen 
number, which, as we will see in the next sub-section, is 
not desirable. 

Just as was found for the Ma number, it is relatively 
speaking more expensive computationally to simulate for 
low Re numbers because the flow velocity must be kept 
low, which means longer simulation times are necessary 
to reach time-scales where the suspended particles or 
fluid flows have moved a significant distance. We there- 
fore compromise and normally keep Re ^ 0.2, which is 
similar to the choice made for LB simulations 0, [t^ ■ 
For many situations related to the flow of colloids, this 
should be sufficiently stringent. 



For large Knudsen numbers Kn ^ 10 continuum Navier- 
Stokes equations completely break down, but even for 
much smaller Knudsen numbers, significant rarefaction 
effects are seen. For example, for flow in a pipe, where a 
in Eq. H34|l is taken to be the pipe radius, important non- 
continuum effects are seen when Kn 0.1. In fact it is 
exactly for these conditions of modest Kn numbers, when 
corrections to the Navier Stokes become noticeable, that 
DSMC approaches are often used 

|3l|33- SRD, which 
is related to DSMC, may also be expected to work well 
for such flows. It could flnd important applications in 
microfluidics and other micromechanical devices where 
Kn effects are expected to play a role js^, • 

Colloidal dispersions normally have very small Kn 
numbers because the mean-free path of most liquid sol- 
vents is very small. For water at standard temperature 
and pressure Afrec ~ 3 A. Just as found for the other di- 
mensionless numbers, coarse-graining typically leads to 
larger Kn numbers because of the increase of the mean- 
free path. Making the Kn number smaller also typically 
increases the computational cost because a larger number 
of collisions need to be calculated. In our simulations, we 
keep Kn< 0.05 for spheres. This rough criterion is based 
on the observation that for small Kn numbers the fric- 
tion coefficient on a sphere is expected to be decreased 
by a factor \ — a Kn, where a is a material dependent 
constant of order 1 [63|, so that we expect Kn number 
effects to be of the same order as other coarse-graining 
errors. There are two ways to achieve small Kn numbers: 
one is by increasing (Tc//ao: the other is by decreasing A. 
The second condition is commensurate with a large Sc 
number or a small Re number. 

In ref. , the Kn numbers for colloids were Kn w 0.5 
and Kn « 0.8, depending on their colloid-solvent cou- 
pling. Such large Kn numbers should have a signiflcant 
effect on particle friction coefficients, and this may ex- 
plain why these authors find that volume fraction has 
less of an effect on the sedimentation velocity than we 
do Isl. 



E. Peclet number 

The Peclet number Pe measures the relative strength of 
convective transport to diffusive transport. For example, 
for a colloid of radius a, travelling at an average velocity 
Vs, the Pe number is defined as: 



D. Knudsen number 

The Knudsen number Kn measures rarefaction effects, 
and can be written as 



Kn 



Afr( 



(34) 



where a is a characteristic length-scale of the fluid flow, 
and Afroc is the mean-free path of the fluid or gas. For a 
colloid in SRD one could take a « Ucf and Afree ~ Aoq. 



Pe 



Vsa 



col 



(35) 



where I?coi is the colloid diffusion coefficient. Just as for 
the Re number, the Pe number can be interpreted as a 
ratio of a diffusive to a convective time-scale, but now the 
former time-scale is not for the diffusion of momentum 
but rather it is given by the colloid diffusion time 



TD 



-Dcol 



(36) 
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which measures how long it takes for a colloid to diffuse 
over a distance CTc/. We again take a ~ CTc/ so that, using 
Eq. H30|l . the Pe number can be written as: 



ts 



(37) 



If Pe ^ 1 then the colloid moves convectively over a dis- 
tance much larger than its radius CTc/ in the time td that 
it diffuses over that same distance. Brownian fluctuations 
are expected to be less important in this regime. For Pe 
<C 1, on the other hand, the opposite is the case, and 
the main transport mechanism is diffusive (note that on 
long enough time-scales {t > td /Pe^) convection will al- 
ways eventually "outrun" diffusion ) . It is sometimes 
thought that for low Pe numbers hydrodynamic effects 
can be safely ignored, but this is not always true. For ex- 
ample, we found that the reduction of average sedimenta- 
tion velocity with particle volume fraction, famously first 
explained by Batchelor |82l. is independent of Pe number 
down to Pe = 0.1 at least l53. 



1. Dependence of Peclet number on simulation parameters 

The highest Pe number achievable in simulation is lim- 
ited by the constraints on the Ma and Re numbers. For 
example the Ma number sets an upper limit on the maxi- 
mum Pe number by limiting Vs- From Eqs. H29|l and (|35|l . 
it follows that the Peclet number can be re-written in 
terms of the Reynolds number as 



Pe 



Dc 



-Re ~ Gttj 



ao 



Re 



(38) 



where we have approximated -Dcoi ~ ksT / {QiTrjUcf). 
This shows that for a given constraint on the Re num- 
ber, increasing 7 or Ucf increases the range of accessi- 
ble Pe numbers. Similarly, when the kinematic viscosity 
is dominated by the coUisional contribution, decreasing 
the dimensionless mean-free path A will also increase the 
maximum Pe number allowed since Pe"*"^ ~ A~^Re ~ 
A^^Ma. However, these changes increase the computa- 
tional cost of the simulation. 
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FIG. 5: Finite system size scaling of the Stokes friction on 
a sphere of size a^f = 2 in a flow of velocity «s = 0.01 ao/to. 
The correction factor / deflned in Eq. 14211 can be extracted 
from measured friction ^ (drag force/velocity) and estimated 
Enskog friction defined in Eq. 14011 . It compares well to 
theory, Eq. CT . 



solved due to Kn number and discretization effects. De- 
creasing the box-size L falls foul of the long-ranged na- 
ture of the HI, and therefore, just as for Coulomb inter- 
actions finite size effects such as those induced by 
periodic images must be treated with special care. 

Increasing the flow velocity may also be desirable since 
more Stokes times ts = Ucf/vs can be achieved for the 
same number of SRD collision steps, thus increasing com- 
putational efficiency. The Ma number gives one con- 
straint on Vg, but usually the more stringent constraint 
comes from keeping the Re number low to prevent un- 
wanted inertial effects. 



A. Finite-size effects 



1. Finite-size correction to the friction 



V. FINITE-SIZE, DISCRETIZATION, AND 
INERTIAL EFFECTS 

The cost of an SRD simulation scales almost linearly 
with the number of fluid particles Nf in the system, and 
this contribution is usually much larger than the cost of 
including the colloidal degrees of freedom. To optimize 
the efficiency of a simulation, one would therefore like 
to keep Nf as small as possible. This objective can be 
achieved by keeping CTc/ /iq and the box-size L small. Un- 
fortunately both these choices are constrained by the er- 
rors they introduce. Reducing u^f ja^ means that short- 
ranged hydrodynamic fields become less accurately re- 



The friction coefficient ^ can be extracted from the 
Stokes drag on a fixed colloid in fluid flow: 

Fd -Cvoo = -47r77aVoo (39) 

where Vqo is the flow fleld at large distances. The pre- 
factor 4 comes from using slip boundary conditions; it 
would be 6 for stick-boundary conditions, as verified 
in 58] . In principle, since 77 is accurately known from 
theory for SRD, this expression can be used to extract 
the hydrodynamic radius a, which is not necessarily the 
same as tJc/ , from a simulation, as we did in [sj . 

The hydrodynamic radius a can also be directly cal- 
culated from theory. To derive this, it is important to 
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recognize that there are two sources of friction |49j . The 
first comes from the local Brownian colhsions with the 
small particles, and can be calculated by a simplified En- 
skog/Boltzmann type kinetic theory |8^ : 



8 (2iTkBTMcm 



E = n 



M 



1/2 



(40) 



which is here adapted for slip boundary conditions. A re- 
lated expression for stick boundary conditions, including 
that for rotational frictions, is described in [53, where it 
was shown that the short-time exponential decays of the 
linear and angular velocity auto-correlation functions are 
quantitatively described by Enskog theory. 

The second contribution to the friction, ^5, comes from 
integrating the Stokes solution to the hydrodynamic field 
over the surface of the particle, defined here as r = cTc/. 
These two contributions to the friction should be added 
in parallel to obtain the total friction ,^83^ .8^ (see also 
Appendix B): 



11 1 



(41) 



In contrast to the Enskog friction, which is local, we 
expect substantial box-size effects on the Stokes friction 
since it depends on long-ranged hydrodynamic effects. 
These can be expressed in terms of a correction factor 
f{<7cf/L) that should go to 1 for very large systems : 



(42) 



To measure this correction factor we plot, in FiglSJ the 
form A-K-qacf (l/'f — I/Cb) for various system box sizes for 
which we have measured ^ from the simulation, and es- 
timated from Eq. H4U|I . As expected, the correction 
factor tends to 1 for smaller acf /L. More detailed calcu- 
lations [HE US J taking into account the effect of periodic 
boundaries, suggests that to lowest order in Ucf / L the 
correction factor should scale as 



f{a,f/L) « 1 + 



2.837^. 

L 



(43) 



Indeed, a least squares fit of the data to this form gives 
a slope of c « 2.9, close to the theoretical value and in 
agreement with similar Lattice Boltzmann simulations of 
a single colloidal sphere . 

With these ingredients in hand, we can calculate the 
theoretical expected friction from Eqns. (|in|l - We 
know from previous work that the Enskog contribution at 
short times and the hydrodynamic contribution at long 
times quantitatively reproduce the rotational and trans- 
lational velocity auto-correlation functions (VACF) for 
fc//ao = 2 — 5 |58j| . and so we expect that the friction 
coefficients should also be accurately described by these 
theories. We test this further in Fig. |H1 for a number of 
different values of acf/iQ, and find excellent agreement 
with theory for Re < 1 and acf/ao > 2. For acf/ao = 1 
and below, on the other hand, we find deviations from 
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FIG. 6: Drag force Fd divided by i-Krjv for various colloid sizes 
and Reynolds numbers (in all cases the box size L = 16(Tc/). 
For small Reynolds numbers this ratio converges to the effec- 
tive hydrodynamic radius a. Dashed lines are theoretical pre- 
dictions from combining Stokes and Enskog frictions in II4H . 



the theory. These are most likely due to Kn number and 
discretization effects, to be discussed in the next sub- 
section. 

For the smallest spheres the mass ratio Mc/mf may 
also change the measured friction if instead of fixing the 
sphere we were to let it move freely. For CTc/ = o-o, 
MrJrri f ~ 21, which is small enough to have a significant 
effect (83 ■ For larger spheres, this is not expected to be 
a problem. For example Mc/mf w 168 for acf = 2ao and 
Mc/mf « 1340 for acf = 4ao. 



2. Effective hydrodynamic radius 

From the calculated or measured friction we can also 
obtain the effective hydrodynamic radius floff, which is 
continuum concept, from a comparison of the microscopic 
frictions from Eq. (glj with Eq. (j^ : 



fleff 



^^ri {^s + ^e) 



(44) 



(for stick-boundaries the 47r should be replaced by Gtt). 
We find acff(a-c/ = 6) = 6.11ao, acsic^cf = 4) = 3.78ao, 
and aofr((Tc/ — 2) = 1.55ao. The effective hydrodynamic 
radius is increased by the finite-size effects, but lowered 
by the Enskog contribution which is added in parallel. 
Because this latter contribution is relatively more impor- 
tant for small colloids, a^s < acf for smaller acf- For 
acf = 6ao the Enskog contribution is smaller than the 
finite-size effects, and so the effective hydrodynamic ra- 
dius is larger than acf- Obviously this latter effect de- 
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pends on the box-size. In an infinite box < Ccf for all 
values of <7c/, due to the Enskog contribution. Note that 
it is the effective hydrodynamic radius Oefi which sets the 
long-ranged hydrodynamic fields, as we will see below. 

3. Turning off long-ranged hydrodynamics 

Hydrodynamic forces can be turned off in SRD by reg- 
ularly randomizing the absolute fluid particle velocities 
(with for example a naive Langevin thermostat). For 
particles embedded in an SRD solvent through partic- 
ipation in the collision step, this trick can be used to 
compare the effects of hydrodynamics to a purely Brow- 
nian simulation|4d|. The Yeomans group has successfully 
applied this idea in a study of polymer collapse and pro- 
tein folding [43. 

However for the case of the colloids embedded through 
direct solvent collisions, turning off the hydrodynamic 
forces by randomizing the velocities greatly enhances the 
friction because, as is clear from Eqns. H4UI) - (|42l) . the 
two contributions add in parallel. Without long-ranged 
hydrodynamics, the friction would be entirely dominated 
by the Enskog contribution (|40|) which scales with a^p 
and can be much larger than the hydrodynamic contri- 
bution which scales as CTc/. Another way of stating this 
would be: By locally conserving momentum, SRD al- 
lows the development of long-ranged hydrodynamic fluid 
velocity correlations that greatly reduce the friction felt 
by a larger colloidal particle compared to the friction it 
would feel from a purely random Brownian heat-bath at 
the same temperature and number density = 7/^0. 

B. Discretization effects on flow- field and friction 
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FIG. 7: Magnified difference fields 10[v(r) — vst(r)]/-i;oo for 
different colloidal sizes (the axes are scaled by (Tc/). In all 
cases Re « 0.1, and box size L = IGctc/. 



with theory for these smallest sphere sizes. We also note 
that using a = Ucf instead of the more accurate value 
of the effective hydrodynamic radius a calculated from 
theory results in significantly larger deviations between 
measured and theoretical flow-fields. This independently 
confirms the values of the effective hydrodynamic radius. 



For pure Stokes flow (Re = 0), the velocity around a 
flxed slip-boundary sphere of (effective) hydrodynamic 
radius a can be exactly calculated: 

vst (r) = (1 - |-) - • f f |- (45) 

where Vqo is the velocity field far away from the sphere, 
and r the vector pointing from the center of the sphere to 
a position inside the fluid, with corresponding unit vector 
f — r/r. 

In Fig. El we plot the difference between the measured 
field and the theoretical expected field H45|) for four dif- 
ferent colloid radius to cell size ratios (Tc//ao, using the 
values of the effective hydrodynamic radius a calculated 
from combining Eq. (|41|l and Eq. As expected, the 
field is more accurately reproduced as the colloid radius 
becomes larger with respect to the cell size ag. This 
improvement arises because SRD discretization and hy- 
drodynamic Knudsen number effects become smaller. In 
Fig. we observe quite large deviations in the hydrody- 
namic field for (Jcf < 1. These discretization effects may 
explain why the measured frictions in Fig.|Sldo not agree 



The increased accuracy from using larger a^f /ao comes 
at a sharp increase in computational cost. To make sure 
that the finite size effects in each simulation of Fig. [7| 
are approximately the same, the box-size was scaled as 
L/ucf- This means that doubling the colloid size leads 
to an eightfold increase in the number of fluid particles. 
Moreover, the maximum velocity of the fluid must go 
down linearly in colloid size to keep the Re number, de- 
fined in Eq. H29|) . constant. If in addition we keep the 
number of Stokes times fixed, meaning that the fiuid fiows 
a certain multiple of the colloid radius or the box size, 
then a larger particle also means the fluid needs to flow 
over a proportionally longer total distance. The over- 
all computational costs for this calculation then scales at 
least as {a■cf/o■o)^^ which is quite steep. For that reason, 
we advocate using smaller colloids wherever possible. 

In most of our simulations we choose tXc/ = 2, which 
leads to a small relative error in the full velocity field and 
for which we can fully explain the observed friction (as 
seen in the previous subsection). This size is similar to 
what is commonly used in LB 0, [t^ . 
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FIG. 8: (a) Fluid flow-field around a fixed colloid with CTc/ = 2 
and L = 64ao for Re = 0.08. (b)-(d) Difference field 10[v(r) - 
vst(r)]/iioo for Re = 0.08, 0.4, and 2, respectively. 



VI. THE HIERARCHY OF TIME SCALES 
A. Colloidal time-scales 

Many different time-scales govern the physics of a col- 
loid of mass M embedded in a solvent 1] . The most im- 
portant ones are summarized in table ITTll and discussed 
in more detail below. 



1. Fluid time-scales 

The shortest of these is the solvent collision time Tcoi 
over which fluid molecules interact with each other. For 
a typical molecular fluid, Tcoi is on the order of a few tens 
of fs, and any MD scheme for a molecular fluid must use a 
discretization time-step Imd ^ TcoI to properly integrate 
the equations of motion. 

The next time-scale up is the solvent relaxation time 
Tf, which measures how fast the solvent VACF decays. 
For a typical molecular liquid, Tf « 10^^^ — 10^^'^ s. 
(For water, at room temperature, for example, it is on 
the order of 50 femtoseconds) 



C. Inertial effects on flow-field and friction 



One of the challenges in SRD is to keep the Re number 
down. It is virtually impossible to reach the extremely 
low Re number (Stokes) regime of realistic colloids, but 
that is not necessary either. In Fig.Elwe see that the fric- 
tion only begins to noticeably vary from the Stokes limit 
at Re « 1 (at least on a logarithmic scale). We expect 
the flow-field itself to be more sensitive to finite Re num- 
ber effects. To study this directly, we examine, in Fig.|Sl 
the flow-field around a fixed colloid of size Cc/ = 2ao in 
a box of L = 32ao for for different values of Re. The dif- 
ferences with the Stokes flow field of Eq. (|45|l are shown 
in Fig.|H|;b)-(d) for Re = 0.08, 0.4, and 2 respectively. In 
all cases we used a hydrodynamic radius of a = 1.55 for 
the theoretical comparison, as explained in the previous 
sub-section. The lengths of the vectors in Figs, (b)-(d) 
are multiplied by 10 for clarity. We observe that the rel- 
ative errors increase with Re. They are on the order of 
5% for Re=0.08, and increase to something on the order 
of 40% for Re=2. This is exactly what is expected, of 
course, as we are moving away from the Stokes regime 
with which these flow lines are being compared. Since 
we also expect effects on the order of a few % from Kn 
number. Ma number and various flnite size effects, we 
argue that keeping Re< 0.2 should be good enough for 
most of the applications we have in mind. 



2. Hydrodynamic time-scales for colloids 

Hydrodynamic interactions propagate by momentum 
(vorticity) diffusion, and also by sound. The sonic time 

tcs = -. (46) 

Cs 

it takes a sound wave to travel the radius of a colloid is 
typically very small, on the order of 1 ns for a colloid of 
radius a = l/xm. For that reason, sound effects are of- 
ten ignored for colloidal suspensions under the assump- 
tion that they will have dissipated so quickly that they 
have no noticeably influence on the dynamics. However, 
some experiments and theory [s^ do flnd effects from 
sound waves on colloidal hydrodynamics, so that this is- 
sue is not completely settled yet. 

The kinematic time, t^, defined in Eq. H31|l (and Ta- 
ble ^Oj as the time for momentum (vorticity) to diffuse 
over the radius of a colloid, is particularly important for 
hydrodynamics. It sets the time-scale over which hy- 
drodynamic interactions develop. For a colloid of radius 
Ifim, T„ K, 10^^ s, which is much faster the colloidal 
diffusion time td. 

When studying problems with a finite fiow velocity, 
another hydrodynamic time-scale emerges. The Stokes 
time ts, defined in Eq. H30|l (and Table lllHl as the time for 
a colloid to advect over its own radius, can be related to 
the kinematic time by the relation = Re tg. Because 
colloidal particles are in the Stokes regime where Re <C 1, 
we find Ti, '^ts- 

Simulation and analytical methods based on the Os- 
een tensor and its generalizations, e.g. 0,0 j implicitly 
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TABLE III: Time-scales relevant for colloidal suspensions 



Solvent time-scales 

Solvent collision time over which solvent molecules interact 

Solvent relaxation time over which solvent velocity correlations decay 



Tcol ~ 10 s 

r/ ~ 10 — 10 s 



Hydrodynamic time-scales 

Sonic time over which sound propagates one colloidal radius 

Kinematic time over which momentum (vorticity) diffuses one colloidal radius 
Stokes time over which a colloid convects over its own radius 



a 

tcs — 

Cs 



a 

V 



ts 



a 

Vs 



TP 

Pe 



Eq. ignj 

Eq. 

Eq. loni 



Brownian time-scales 

Fokker-Planck time over which force-force correlations decay 

Enskog relaxation time over which short-time colloid velocity correlations decay 



tfp 



TE 



Brownian relaxation time over which colloid velocity correlations decay in Langevin Eq. tb = 



Colloid diffusion time over which a colloid diffuses over its radius 



Eq. m 
Eq. igTJ 
Eq. lOni 



Ordering of time-scales for colloidal particles 



Tcol < Tf,TFP < TEjtcs < TB < Ty < TD,ts 



assume that the hydrodynamic interactions develop in- 
stantaneously, i.e. that « 0. For the low Re number 
regime of colloidal dispersions this is indeed a good ap- 
proximation. Of course it must be kept in mind that t^^ 
is a diffusive time-scale, so that the distance over which 
it propagates grows as \/t. Thus the approximation of 
instantaneous hydrodynamics must be interpreted with 
some care for effects on larger length-scales. 

3. Brownian time-scales for colloids 

The fastest time-scale relevant to the colloidal Brown- 
ian motion is the Fokker-Planck time Tpp defined in ij as 
the time over which the force-force correlation function 
decays. It is related to rj, the time over which the fluid 
looses memory of its velocity because the forces on the 
colloid are caused by collisions with the fluid particles, 
velocity differences de-correlate on a time scale of order 

The next time-scale in this series is the Brownian time: 



where ^5 = Q^rja is the Stokes friction for stick bound- 
ary conditions. It is often claimed that this measures the 



time for a colloid to lose memory of its velocity. How- 
ever, as discussed in Appendix^ this picture, based on 
the Langevin equation, is in fact incorrect for colloids. 
Nevertheless, tb has the advantage that it can easily be 
calculated and so we will use it as a crude upper bound 
on the colloid velocity de-correlation time. 

Perhaps the most important time-scale relevant for 
Brownian motion is the diffusion time td, described in 
Eq. H3t)|) as the time for a colloidal particle to diffuse 
over its radius. For a colloid of radius 1/xm, r^i w 5 s, 
and even though td ~ a'^, it remains much larger than 
most microscopic time-scales in the mesoscopic colloidal 
regime. 



B. Time-scales for coarse-grained simulation 

In the previous subsection we saw that the relevant 
time-scales for a single colloid in a solvent can span as 
many as 15 orders of magnitude. Clearly it would be 
impossible to bridge all the time-scales of a physical col- 
loidal system - from the molecular TcoI to the mesoscopic 
T£) - in a single simulation. Thankfully, it is not neces- 
sary to exactly reproduce each of the different time-scales 
in order to achieve a correct coarse-graining of colloidal 
dynamics. As long as they are clearly separated, the cor- 
rect physics should still emerge. 
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In SRD the physical time-scale Tcoi is coarse-grained 
out, and the effect of the collisions calculated in an aver- 
age way every time-step A^c- The time-scale Tf on which 
the velocity correlations decay can be quite easily calcu- 
lated from a random-collision approximation. Following 
H^: r/ w -Aio/ln[l-|(l-cosa)(l-l/7)]. For our pa- 
rameters, a = ^TT, 7 = 5, we find Tf « 0.76Atc — 0.076io- 
However, it should be kept in mind that for small A the 
exponential decay with t/r/ turns over to a slower alge- 
braic decay at larger t [Sg- The amplitude of this "tail" 
is, however, quite small. 



FIG. 9: (Color online) Schematic depiction of our strategy 
for coarse-graining across the hierarchy of time-scales for a 
colloid (here the example taken is for a colloid of radius 1 /im 
in H2O). We first identify the relevant time-scales, and then 
telescope them down to a hierarchy which is compacted to 
maximise simulation efficiency, but sufficiently separated to 
correctly resolve the underlying physical behaviour. 



Since we take the view that the "solvent" particles are 
really a Navier Stokes solver with noise, what is needed 
to reproduce Brownian behavior is first of all that the col- 
loid experiences random kicks from the solvent on a short 
enough time-scale. By dramatically reducing the num- 
ber of "molecules" in the solvent, the number of kicks 
per unit of time is similarly reduced. Here we identify 
the Fokker-Planck time Tpp as the time-scale on which 
the colloid experiences random Brownian motion. For 
Brownian motion it doesn't really matter how the kicks 
are produced, they could be completely uncorrelated, but 
since we also require that the solvent transports momen- 
tum, solvent particles must have some kind of correla- 
tion, which in our case is represented by the time Tf. 
Other colloid relaxation times should be much longer 
than this. We therefore require first of all that Tpp tb 
and Tf ^ tb- 

Secondly, the colloid's VACF should decay to zero well 
before it has diffused or convected over its own radius. 
A proper separation of time-scales in a coarse-graining 
scheme then requires that tb <C tb as well as tb <^ ts 
for systems where convection is important. 

Finally, for the correct hydrodynamics to emerge we 
require that Tpp,Tf <C ^ tb- But this separation no 
longer needs to be over many orders of magnitude. As ar- 
gued in [s^ l , one order of magnitude separation between 
time-scales should be sufficient for most applications. We 
illustrate this approach in Fig. |^ 

In the next few subsections we discuss in more quanti- 
tative detail how our coarse-graining strategy telescopes 
down the hierarchy of time-scales in order to maximize 



2. Hydrodynamic time-scales for simulation 

For our choice of units Cs = -\/5/3 ao/to, so that the 
sonic time of Eq. H4ti|) reduces to 



tc 



0.775' 



ao ' 



(48) 



which is independent of A or 7. 

In the limit of small A, the ratio of the kinematic time 
to Ate can be simplified to the following form: 



18- 



(49) 



so that the condition Tf,Tpp <C Ti, is very easy to fulfil. 
Furthermore, under the same approximations, the ratio 
T^/tcs ~ 23(Tc/A so that for A too small the kinematic 
time becomes faster than the sonic time. For the simu- 
lation parameters used in |52| | (and in Fig. 110(1 . we find 

Tv j Tcs ~ 5 . 



3. Brownian time-scales for simulation 

We expect the Fokker Planck time Tpp to scale as Ate, 
since this is roughly equivalent to the time Tf over which 
the fluid velocities will have randomized. In Fig. EH we 
plot the force-force correlation function 



Cpp{t) 



{F{t)Fm 
in ■ 



(50) 



Its short time behavior is dominated by the random 
forces, and the initial decay time gives a good estimate of 
Tpp. As can be seen in the inset of Fig. ^1 Tpp « 0.09io, 
which is indeed on the order of At^ or Tf. 

If we make the reasonable approximation that cTcc ~ 
2crcf, then the ratio of the Brownian time tb to the kine- 
matic time simplifies to: 

TB_ _ 2pc 



(51) 
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FIG. 10: (Color online) Normalized VACF for a single col- 
loid with (Tc/ — 2ao and L — 32ao: measured (solid line), 
Enskog short time prediction (dot-dashed line), and Brown- 
ian approximation (dashed line). The normalized force- force 
auto-correlation Cppit) (dotted line) decays on a much faster 
time scale. The inset shows a magnification of the short-time 
regime. Note that the colloid exhibits ballistic motion on a 
time-scale ~ r/. The time-scales tfp ~ 0.09, tb ~ 2.5 and 
td ~ 200 are clearly separated by at least an order of magni- 
tude, as required. 



SO that the ratio of tb to is the same as for a real 
physical system. Since buoyancy requirements mean that 
normally pc ~ p/, the time-scales are ordered as tb < t^- 
Moreover, since the ratio r^/Atc ^ 1, independently of 
7 or even A (as long as A <C 1), the condition that Tpp <C 
Tb is also not hard to fulfil in an SRD simulation. 

These time-scales are illustrated in Fig. ^| where we 
plot the normalized time correlation function: 



'-^VVW - : 



(52) 



with V the Cartesian velocity coordinate of a single col- 
loid. After the initial non-Markovian quadratic decay, 
the subsequent short-time exponential decay is not given 
by Tb but instead by the Enskog time 



TE = 



Me 



(53) 



where is the Enskog friction which can be calculated 
from kinetic theory see Eq. (|^ . and generally 

te < Tb- The physical origins of this behavior are de- 
scribed in more detail in Appendix IbI 

The colloid diffusion coefficient is directly related to 
the friction by the Stokes-Einstein relation: 



kpT 



(54) 



FIG. 11: Time dependent self-diffusion coefficient D{t) for 
different volume fractions (p, achieved by varying the number 
of colloids in a box of fixed volume. The inset shows the 
self-diffusion coefficient {D — \imt ooD{t) as a function of 
volume fraction, together with an analytical prediction D = 
L>o(l - 1.1795(^) ^ valid in the low-density limit. 



If we assume that A ^ 1 and, for simplicity, that ^ ~ Csi 
i.e. we ignore the Enskog contribution, then the diffusion 
time scales as: 



TD 



D. 



col 



OcS_ 



7 
4A2 



ao 



TB, 

(55) 

so that Tb ^ tu is not hard to fulfil. 

It is also instructive to examine the ratio of the diffu- 
sion time Tb, to the kinematic time r^: 



Tu 



A^ V «o 



(56) 



In general we advocate keeping A small to increase the Sc 
number Sc = i^/Df, and since another obvious constraint 
is Dcoi Df, there is not too much difhculty achieving 
the desired separation of time-scales t^ <ti td- 

As a concrete example of how these time-scales are 
separated in a simulation, consider the parameters used 
in Fig.Hnifor which we find t/ = 0.076io, Tpp = 0.09to 
Tb — 2.5io, — 8tQ and Tp = 200io. But more generally, 
what the analysis of this section shows is that obtaining 
the correct hierarchy of time-scales: 



Tf,Tpp < Te,Tb < T^ <^TD,ts 



(57) 



is virtually guaranteed once the conditions on dimension- 
less numbers, detailed in Section Hvl are fulfilled. 
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FIG. 12: Mean-square-displacement of a ctc/ = 2ao colloid 
(solid line). The dotted line is a fit to the long time limit. The 
inset shows a magnification of the short-time regime, high- 
lighting the initial ballistic regime (dashed line is (V^(0))^f^), 
which rapidly turns over to the diffusive regime. 



C. Measurements of diffusion 

To further test the hybrid MD-SRD coarse-graining 
scheme, we plot in Fig. 1111 the time-dependent diffusion 
coefficient D(t\ defined as: 

D{i)^ fdt'{V{t')vm (58) 
Jo 

for a number of different volume fractions (j) = 
/Vna'^j. Note that the convergence slows with in- 
creasing volume fraction (p. The infinite time integral 
gives the diffusion coefficient, i.e. Z^coi = fimt^txj £>(<). 
In the limit of low densities the diffusion coefficient has 
been predicted to take the form = Dq{1 — c0) with 

1.1795 for slip boundary spheres |90|, and this provides 
a good fit at the lower volume fractions. 

In Figure IT^ we plot the mean-square displacement of 
a Cartesian component of the colloidal particle position. 
As highlighted in the inset, at short times there is an 
initial ballistic regime due to motion at an average mean 
square velocity (^(0)^) = ksT/M, resulting m a mean- 
square displacement {{X{t)-X{Q))'^) = {kBTlM)t^. At 
times t ^ Ti,, the motion is clearly diffusive, as expected, 
with a linear dependence on t and a slope of 2Dcoi, where 
Dcoi is given by the infinite time limit of Eq. (|58|1 . On 
the scale of the plot the asymptotic regime is reached well 
before the time td over which the particle has diffused, 
on average, over its radius. In Appendix IbI the behavior 
of D(t) and the related mean-square displacement are 
discussed in more detail for time-scales t ^ to- 



VII. MAPPING BETWEEN A PHYSICAL AND 
COARSE-GRAINED SYSTEMS: NO FREE 
LUNCH? 

The ultimate goal of any coarse-graining procedure is 
to correctly describe physical phenomena in the natural 
world. As we have argued in this paper, it is impossible 
to bridge, in a single simulation, all the time and length- 
scales relevant to a colloid suspended in a solvent. Com- 
promises must be made. Nevertheless, a fruitful mapping 
from a coarse-grained simulation to a physical system is 
possible, and greatly facilitated by expressing the phys- 
ical properties of interest in dimensionless terms. The 
best way to illustrate this is with some examples. 



A. Example 1: Mapping to diffusive and Stokes 
time-scales 

For many dynamic phenomena in colloidal dispersions, 
the most important time-scale is the diffusion time to- 
To map a coarse-grained simulation onto a real physi- 
cal system one could therefore equate the diffusion time 
and the colloid radius a of the simulation to that of the 
real physical system. For example, take the simulation 
parameters from table IVll in Appendix lO for a = 2ao, 
and compare these to the properties of colloids of radius 
a = 1/im and a = 10 nm from tabled For the larger 
colloids. To — 5s, and equating this to the simulation to 
means that oq = 0.5/.tm and to = 0.02 s. Performing the 
same exercise for the smaller colloid, where r^, = 5 x 10~^ 
s, results in ao = 5 nm and to — 2x 10~* s. The first thing 
that becomes obvious from this exercise is that "physi- 
cal" time is set here not by the coarse-grained simulation, 
but by the particular system that one is comparing to. 
In other words, many different physical systems could be 
mapped to the same simulation. 

If the system is also subjected to flow, then a second 
time-scale emerges, the Stokes time ts = Pe to- As long 
as the physical Pe number is achievable without compro- 
mising the simulation quality, then this time-scale can si- 
multaneously be set equal to that of the physical system. 
Behavior that depends primarily on these two time-scales 
should then be correctly rendered by the SRD hybrid MD 
scheme described in this paper. 



B. Example 2: Mapping to kinematic time-scales 

By fixing To or ts, as in example 1, other time-scales 
are not correctly reproduced. For the large and small col- 
loid systems described above, we would find t,^ — 0.17 s 
and 1.7 x 10~^ s respectively, which contrasts with the 
physical values of — 10~^ s and 10^^" s shown in Ta- 
ble In other words, fixing to results in values of t^ 
that are too large by orders of magnitude. But of course 
these much larger values of t^, are by design, because 
the SRD simulation is optimized by compacting the very 
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large physical hierarchy of time-scales. If we are inter- 
ested in processes dominated by td, having different 
doesn't matter, so long as <^td- 

On the other hand, if we want to describe processes 
that occur on much smaller time-scales, for example long- 
time tails in colloidal VACFs, then simulation times Ti, 
could be mapped onto the physical instead. For the 
same a — 2ao simulation parameters used above, the 
a = 1/im system now maps onto ag — 0.5/im, and 
io = 1.3 X 10^^ s with of course a strongly underesti- 
mated diffusion time: tq — i x 10~^ s instead of 5 s. 
Similarly for the a = 10 nm system, the mapping is to 
flo = 5 nm, to — 1-3 x 10~^^ s and td = 3 x 10~® s instead 
of 5 X 10^^ s. For many processes on time-scales t < td, 
this underestimate of tjj doesn't really matter. It is only 
when time-scales are mixed that extra care must be em- 
ployed in interpreting the simulation results. And, if the 
processes can be described in terms of different dinien- 
sionless times, then it should normally still be possible 
to disentangle the simulation results and correctly map 
onto a real physical system. 

Another lesson from these examples of mapping to dif- 
ferent time-scales comes from comparing the kinematic 
viscosity of a physical system, which say takes the value 
1/ — 10^/im^/s for water, to its value in SRD simula- 
tions. If for the a = 2ao SRD system described above 
we map time onto td, as in Example 1, then for the 
a = 1/im system v = 6/im^/s, while for the a = 10 nm 
system v — Qx 10~^/im^/s. As expected, both are much 
smaller than the true value in water. If we instead map 
the times to t^^ then of course the kinematic viscosity 
takes the same value as for the physical system (since 
we also set the length-scales equal to the physical length- 
scales). This apparent ambiguity in defining a property 
like the kinematic viscosity is inherent in many coarse- 
graining schemes. We take the view that the precise value 
of V is not important, the only thing that matters is that 
it is large enough to ensure a proper separation of times- 
scales and the correct regime of hydrodynamic numbers. 



C. Example 3: Mapping mass and temperature 

In the examples above, we have only discussed setting 
lengths and times. This leaves either the mass m/ or tem- 
perature (thermal energy) fc^T still to be set. Because for 
many dynamic properties the absolute mass is effectively 
an irrelevant variable (although relative masses may not 
be) there is some freedom to choose. One possibility 
would be to set Mc, the mass of the colloid, equal to 
the physical mass. For an a = 1/xm neutrally buoyant 
colloid this sets ruf — 2.5 x 10"^'' g, (equal to about 
8 X 10^ water molecules) while for an a = 10 nm colloid 
we find m/ = 2.5 X 10^^" g (equal to about 800 water 
molecules) . By construction this procedure produces the 
correct physical pf. If we fix time with r^, this will there- 
fore also generate the correct shear viscosity rj. If instead 
we fix T£), then rj will be much smaller than the physical 



value. For either choice of time-scales, other properties, 
such as the fluid self-diffusion constant or the number 
density, will have different values from the underlying 
physical system. 

Setting the mass in this way means that the unit of 
thermal energy fc^To = m/aQ^Q^ would also be very 
large, leading to the appearance of very low tempera- 
tures. However, because temperature only determines 
behavior through the way it scales potentials, it is quite 
easy to simply re-normalize the effective energy scales 
and obtain the same behavior as for the physical sys- 
tem. Alternatively one could set the temperature equal 
to that of a physical system, but that would mean that 
the masses would again differ significantly from the real 
physical system. At any rate, for the same simulation 
(say of sedimentation at a given value of Pe) the values 
of mass or temperature depend on the physical system 
one compares to. This apparent ambiguity (or flexibil- 
ity) simply reflects the fact that the dominant effects of 
temperature and mass come into play as relative ratios 
and not as absolute values. 



D. Example 4: Mapping to attractive potentials: 
problems with length-scales? 

So far we have only treated one length-scale in the 
problem, the radius a. Implicitly we have therefore as- 
sumed that the colloid-colloid interaction does not have 
another intrinsic length-scale of its own. For hard-sphere 
colloids this is strictly true, and for steep-repulsions such 
as the WCA form of Eq. (|16|l it is also a fairly good 
assumption that the exact choice of the exponent n in 
this equation does not signiflcantly affect the dynamical 
properties [Blj . Similarly, there is some freedom to set 
the repulsive fluid-colloid potential of Eq. (|17|l in such a 
way as to optimize the simulation parameters, in par- 
ticular AtMD- The physical liquid-colloid interaction 
would obviously have a much more complex form, but 
SRD coarse-grains out such details. A similar argument 
can be made for colloid and fluid interactions with hard 
walls, needed, for example, to study problems with con- 
finement, for which, inter alia, SRD is particularly well 
suited. 

Things become more complicated for colloids with an 
explicit attractive interaction, such as DLVO or a de- 
pletion potential 0. These potentials introduce a new 
length-scale, the range of the attraction. The ratio of 
this length to the hard-core diameter helps determine 
the equilibrium properties of the fluid [23, In a 

physical system this ratio may be quite small. Keep- 
ing the ratio the same in the simulations then leads to 
potentials that are very steep on the scale of acc- The 
MD time-step AImd may need to be very small in or- 
der to properly integrate the MD equations of motion. 
Such a small AImd can make a simulation very ineffi- 
cient, as was found in [s^l who followed this strategy 
and were forced to use AImd/ Ate ~ 455. However, in 
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that case the DLVO potential dominates the behaviour, 
so that SRD was only included to roughly resolve the long 
range hydrodynamics In general, we advocate adapt- 
ing the potential to maximize simulation efficiency, while 
preserving key physical properties such as the topology 
of the phase diagram. For example, the attractive en- 
ergy scale can be set by the dimensionless reduced second 
virial coefficient, which gives an excellent approximation 
for how far one is from the liquid- liquid critical point [92| . 
When the attractive interaction range is less than about 
30% of the colloid hard-core diameter, the colloidal "liq- 
uid" phase becomes metastable with respect to the fluid- 
solid phase-transition. At low colloid packing fraction (/)c 
this leads to so-called "energetic fluid" behavior . To 
simulate a colloidal suspension in this "energetic fluid" 
regime, having a range of less than 30% of a should suf- 
fice. In this way adding attractions does not have to make 
the simulation much less efficient. Moreover the effects 
of these constraints, placed by the colloid-colloid interac- 
tion on AtuD, are to some degree mitigated by the fact 
that it is usually the colloid-fluid interaction which sets 
this time-scale. 

There are still a few final subtleties to mention. First 
of all, using small acf/o-o in the simulations may greatly 
increase efficiency, but it also leads to a relatively larger 
contribution of the Enskog friction to the total friction 
in Eq. H41() . For physical colloids in a molecular solvent 
the exact magnitude of the Enskog friction is still an 
open question, but undoubtedly for larger colloids it is 
almost negligible. Some of the Enskog effects can simply 
be absorbed into an effective hydrodynamic radius a, but 
others must be interpreted with some care. 

Another regime where we might expect to see devia- 
tions beyond a simple re-scaling of time-scales would be 
at very short times. For example, when colloids form a 
crystal, their oscillations can be de-composed into lattice 
phonons. All but the longest wave-length longitudinal 
oscillations are overdamped due to the viscous drag and 
backflow effects of the solvent 93|. In SRD, however, 
some of the very short wavelength oscillations might be 
preserved due to the fact that the motion on those time- 
scales is still ballistic. These may have an impact on the 
interpretation of SRD simulations of microrheology. 

In summary, these examples demonstrate that as long 
as the coarse-grained simulation produces the correct hy- 
drodynamic and Brownian behavior then there is consid- 
erable freedom in assigning real physical values to the pa- 
rameters in the simulation. Attractive interactions may 
introduce a new length-scale, but a judicious choice of a 
model potential should still reproduce the correct phys- 
ical behavior. The same simulation may therefore be 
mapped onto multiple different physical systems. Con- 
versely, for the same physical system, different choices 
can be made for the time-scales that are mapped to phys- 
ical values, but by design, not all time-scales can be si- 
multaneously resolved. 



VIII. CONCLUSION 

Correctly rendering the principal Brownian and hy- 
drodynamic properties of colloids in solution is a chal- 
lenging task for computer simulation. In this article we 
have explored how treating the solvent with SRD, which 
coarse-grains the collisions between solvent particles over 
both time and space, leads to an efficient solution of the 
thermo-hydrodynamic equations of the solvent, in the 
external field provided by the colloids. Although it is 
impossible to simulate the entire range of time-scales en- 
countered in real colloidal suspensions, which can span 
as much as 15 orders of magnitude, we argue that this 
is also not necessary. Instead, by recognizing firstly that 
hydrodynamic effects are governed by a set of dimen- 
sionless numbers, and secondly that the full hierarchy 
of time-scales can be telescoped down to a much more 
manageable range while still being properly physically 
separated, we demonstrate that the key Brownian and 
hydrodynamic properties of a colloidal suspension can 
indeed be captured by our SRD based coarse-graining 
scheme. 

In particular, to simulate in the colloidal regime, the 
modeler should ensure that the dimensionless hydrody- 
namic numbers such as the Mach, Reynolds and Knud- 
sen numbers are small enough I), and the Schmidt 
number is large enough (3> 1). While these numbers are 
constrained by approximate limits, the Peclet number, 
which measures the relative strength of the convective 
transport to the diffusive transport of the colloids, can 
be either larger or smaller than 1, depending on physi- 
cal properties such as the colloidal size and the strength 
of the external field that one wants to study. It turns 
out that if the hydrodynamic numbers above are chosen 
correctly, it is usually not so difficult to also satisfy the 
proper separation of the time-scales. 

We explored how to reach the desired regime of hydro- 
dynamic numbers and time-scales by tuning the simula- 
tion parameters. The two most important parameters are 
the dimensionless mean-free path A which measures what 
fraction of a cell size oq an SRD particle travels between 
collisions performed at every time-step A.tc, and the ratio 
of the colloid radius Uc/ to the SRD cell size ao- The gen- 
eral picture that emerges is that the collision time Ate 
must be chosen so that A is small since this helps keep 
the Schmidt number high, and both the Knudsen and 
Reynolds numbers low. Of course for computational effi- 
ciency, the collision time should not be too small either. 
Similarly although a larger colloid radius means that the 
hydrodynamic fields are more accurately rendered, this 
should be tempered by the fact that the simulation ef- 
ficiency drops off rapidly with increasing colloid radius. 
We find that choosing CTc/ /oo ~ 2, which may seem small, 
already leads to accurate hydrodynamic properties. 

A number of subtleties occur at the fluid-colloid inter- 
face, which may induce spurious depletion interactions 
between the colloids, but these can be avoided by care- 
ful choice of potential parameters for slip boundary con- 
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ditions. Stick boundary conditions can also be imple- 
mented, and we advocate using stochastic bounce-back 
rules to achieve this coarse-graining over the fluid-colloid 
interactions. 

The simplicity of the SRD solvent facilitates the calcu- 
lation of the the short-time behavior of the VACF from 
kinetic theory. We argue that for times shorter than the 
sonic time tcs — a/cg, where hydrodynamic collective 
modes will not yet have fully developed, the decay of the 
VACF is typically much more rapid than the prediction 
of the simple Langevin equation. The ensuing Enskog 
contribution to the total friction should be added in par- 
allel to the microscopic hydrodynamic friction, and from 
the sum an analytic expression for the effective hydrody- 
namic radius aoff can be derived which is consistent with 
independent measurements of the long-ranged hydrody- 
namic fields. 

Although we have shown how to correctly render the 
principal Brownian and hydrodynamic behavior of col- 
loids in solution there is, as always, no such thing as a 
free lunch. The price to be paid is an inevitable con- 
sequence of compressing together the hierarchy of time 
and length-scales: not all the simulation parameters can 
be simultaneously mapped onto a physical system of in- 
terest. However the cost of the "lunch" can be haggled 
down by recognizing that only a few physical parameters 
are usually important. For example, one could map the 
simulation onto real physical time, say through td and 
ts, or alternatively through t^; there is some freedom 
(or ambiguity) in the choice of which time scale to use. 
The correspondence with physical reality becomes more 
complex when different physical time-scales mix, but in 
practice they can often be disentangled when both the 
coarse-grained simulation and the physical system are ex- 
pressed in terms of dimensionless numbers and ratios. In 
other words, there is no substitute for careful physical 
insight which is, as always, priceless. 

A number of lessons can be drawn from this exercise 
that may be relevant for other coarse-graining schemes. 
SRD is closely related to Lattice-Boltzmann and to the 
Lowe- Anderson thermostat, and many of the conclusions 
above should transfer in an obvious way. For dissipative 
particle dynamics, we argue that the physical interpreta- 
tion is facilitated by recognizing that, just as for SRD, the 
DPD particles should not be viewed as "clumps of fluid" 
but rather as a convenient computational tool to solve the 
underlying thermo-hydrodynamic equations. All these 
methods must correctly resolve the time-scale hierarchy, 
and satisfy the relevant hydrodynamic numbers, in order 
to reproduce the right underlying physics. 

Our measurements of the VACF and the discussion 
in Appendix B show explicitly that Brownian Dynam- 
ics simulations, which are based on the simple Langevin 
equation, do not correctly capture either the long or the 
short-time decay of colloidal VACFs. 

A major advantage of SRD is the relative ease with 
which solutes and external boundary conditions (walls) 
are introduced. Boundaries may be hard or soft, with 



either stick or slip conditions. This method is therefore 
particularly promising for simulations in the fields of bio- 
and nanofiuidics, where small meso particles flow in con- 
strained geometries, possibly confined by soft fluctuating 
walls. We are currently exploring these possibilities. 
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APPENDIX A: THERMOSTATING OF SRD 
UNDER FLOW 



When an external field is applied to the fluid (or to ob- 
jects embedded in the fluid), energy is pumped into the 
system and, as a consequence, the average temperature 
will rise. To prevent this from happening, the system 
must be coupled to a thermostat. In order not to influ- 
ence the average flow, thermostating requires a local and 
Galilean invariant definition of temperature. 

We achieve this by relating the instantaneous local 
temperature in a cell to the mean square deviation of the 
fluid particle velocities from the center of mass velocity 
of that cell. To minimize interference of the thermostat 
with the dynamics, we choose to measure the overall tem- 
perature as: 



EE(^»-^"n,c)' (Al) 

cell c iGc 

3(iVccUc-l) (iVcollc> 1) 



iVfrcc - ] 
cell c 







(iVeellc < 1) 



(A2) 



In Eq. (|A2|) , A^ccii c is the instantaneous number of fluid 
particles within cell c. Three degrees of freedom must be 
subtracted for fixing the center of mass velocity of each 
cell (note that this implies that the local temperature is 
not defined in cells containing or 1 particle). During 
the simulations, the temperature Tmoas is measured every 
Ate (this can be done very efficiently because the relative 
velocities are readily available in the collision step rou- 
tine). The thermostat then acts by rescaling all relative 
velocities — Vcm.c by a factor ^/T/Tmcas- This strict 
enforcement of the overall temperature may be relaxed 
by allowing appropriate fluctuations, but in view of the 
very large number of fluid particles (typically 10^) this 
will not have a measurable effect on the dynamics of the 
fluid particles. 
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APPENDIX B: THE L ANGEVIN EQUATION 
AND MEMORY EFFECTS 



Within the Brownian approximation each Cartesian 
component of the coUoid velocity V is described by a 
simple Langevin equation of the form: 

Mc^ - -CV + F^'it) (Bl) 
{F"it)F^{t')) = 2kBT^S{t~t') (B2) 

without any memory effect in the friction coefficient ^, 
or in the random force F^{t). This Langevin equation 
fundamentally arises from the assumption that the force- 
kicks are Markovian: each step is independent of any pre- 
vious behavior. Of course on a very short time-scale this 
is not true, but since Brownian motion is self-similar, the 
Markovian assumptions underlying Eq. (jBip might be ex- 
pected to be accurate on longer time-scales. Inspection 
of Fig. ^1 shows that the force-force correlation func- 
tion (|5U|I does indeed decay on a time-scale Tpp which 
is much faster than other Brownian time-scales in the 
system. 

From the Langevin Eq. (|B1|I it follows that the velocity 
auto-correlation function (VACF) takes the form 

{Vit)Vm'^='^eM-t/rB) (B3) 

where tb = M^/^. The VACF is directly related to the 
diffusion constant through the Green-Kubo relation: 



Dco\ = hm 

t — >oo 



dt{v{t)vm^ 



(B4) 



and Eq. ljB3|) has the attractive property that its integral 
gives the correct Stokes-Einstein form for Dcoi, and that 
its i = limit gives the expected value from equipar- 
tition. It is therefore a popular pedagogical device for 
introducing Brownian motion. Notwithstanding its se- 
ductive simplicity, we will argue that this Langevin ap- 
proach strongly underestimates the initial decay rate of 
the VACF, and badly overestimates its long-time decay 
rate. 



1. VACF at long times 

At longer times, as first shown in MD simulations by 
Alder and Wainwright '93| , the VACF shows an algebraic 
decay that is much slower than the exponential decay of 
Eq. (jB3|) . The difference is due to the breakdown of the 
Markovian assumption for times t 3> Tpp. Memory ef- 
fects arise because the momentum that a colloidal parti- 
cle transfers to the fluid is conserved, with dynamics de- 
scribed by the Navier Stokes equations, resulting in long- 
ranged flow patterns that affect a colloid on much longer 
time-scales than Ea. (jB2|l suggests. Calculations taking 
into account the full time-dependent hydrodynamics were 



performed by Hauge 1951 and Hinch [9g (and others ap- 
parently much earlier |93). For pc = Pf, i.e. a neutrally 
buoyant particle, these expressions can be written as 




— / dx 



a; 2 exp [-xft/r^ )] 



where we have chosen to use an integral form as in 
Under the assumption that the sound effects have dis- 
sipated away, the only hydrodynamic time-scale is the 
kinematic time t^, defined in Eq. 1)3 The VACF must 
therefore scale with i/rj^, independently of colloid size. 
Eq. IjBSp is written in a way to emphasize this scaling. 



By expanding the denominator it is not hard to see 
that at longer times the well-known algebraic tail results: 



lim {V{t)V{0))" 



Mc 9y/n \ t 



12m7 {TTvif^ 
(B6) 

with the same pre- factor that was found from mode- 
coupling theory |99| . 



(B7) 



The integral of Eq. HB5|) is 



which means that all the hydrodynamic contributions to 
the friction f come from {V{t)V{Q))^ . The integral con- 
verges slowly, as {t/T^)^, because of the long-time tail. 
The VACF can easily be related to the mean-square dis- 
placement, and deviations from purely diffusive motion 
are still observable for t ^ r ,^ , and have been measured 
in experiments lol lToollloH . 



Although the last equality in Eq. I|B6(I shows that the 
VACF tail amplitude is independent of colloid size, it 
should be kept in mind that the algebraic decay does not 
set in until t ~ a- /v. Thus it is the relative (not 

absolute) effect of the tail that is the same for all colloids. 
For example, at t — t^, the VACF will have decayed to 
4% of its initial value of ksT/Mc, independently of col- 
loid mass Mc- Even though this tail amplitude may seem 
small, about half the weight of the integral HB7|I that de- 
termines the friction is from t ^ t^. The dominance of 
the non-Markovian hydrodynamic behavior in determin- 
ing the friction for colloids does not appear to hold for 
microscopic fluids, where, for example for LJ liquids near 
the triple point, it is thought that the long-time tail only 
adds a few % to the friction coefficient |l02j . However, 
this is not due to the simple Langevin equation (|Bip be- 
ing more accurate, but rather to the modified short-time 
behavior of the VACF, caused by fluid correlations (a 
different non-Markovian effect). 
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FIG. 13: (Color online) Normalized VACF for a = 2ao 
and stick boundary conditions (from [H^), compared to the 
short-time Enskog result IjBll^ . and the full hydrodynamic 
VACF by itself, and with the VACF lEgJ from sound 

added. We also show the asymptotic long-time hydrodynamic 
tail from mode-coupling theory 



2. VACF at short times 

At short times the hydrodynamic contribution to the 
VACF JbU reduces to 



,i_„,(v(„v(0))» . 



(B8) 



because, within a continuum description, ^ of the energy 
is dissipated as a sound wave at velocity Cg, 



{v{t)vio)Y 



IfceT 
3 M, 



cos 



exp [-{3/2){t/tcs)] ■ 
'#^Vv3sinf^A^ 



(B9) 

where the sonic time t^s = a/cg is defined in Table ITTTl 
The sound wave doesn't contribute to the friction or dif- 
fusion since 



dt{v{t)v{o)y' = 0. 



(BIO) 



For a more extended discussion of the role of sound on 
hydrodynamics see e.g. ref. @. 

Because tcs ^ the sound-wave contribution to the 
VACF decays much faster than the hydrodynamic contri- 
bution {V{t)V{0))^ or even than the Langevin approxi- 
mation {V{t)V{0)) (recall that t^, — for neutrally 
buoyant colloids.) Therefore, in the continuum picture, 
the short-time decay of the VACF is much faster than 



that suggested by the simple Langevin equation. And 
even after the sound wave has decayed, {V{t)V{0)r de- 
cays much more rapidly (faster than simple exponential) 
than Eq. (jB3|) for times t < tb- SRD simulations con- 
firm this more rapid short time decay. For example, in 
Fig 1131 we show a direct comparison of the measured 
VACF to the continuum approach, where {V{t)V{0)) = 
{V{t)V{{)))" + {V{t)V{Q)f\ For short times the decay 
predicted from sound agrees reasonable well, and for long 
times the full hydrodynamic theory quantitatively fits the 
simulation data. 

A closer look at the short-time behavior in our SRD 
simulations, both for slip boundaries as in Fig. ^| and 
for stick-boundaries as in Fig ^| or in js^l , shows that 
a better fit to the short-time exponential decay of the 
VACF is given by 



where the Enskog time (|53|l scales as: 



TE 



0.5 a 



VkeT/r 



0M5tr 



(Bll) 



(B12) 



for heavy (stick boundary) colloids in an SRD fluid. Note 
that the Enskog decay time is very close to the decay 
time |fcs of the sound wave (|B9|) . That the SRD Enskog 
time should scale as tcs is perhaps not too surprising, 
because that is the natural time-scale of the SRD fluid. 
However, the fact that the pre-factors in the exponen- 
tial of Eqs (|B9|I and IjBlip are virtually the same may 
be accidental. Although in a real liquid the increased 
compressibility factor means te and Tcs are lower than 
their values for an ideal gas, the two times change at 
slightly different rates so that we don't expect the same 
equivalence in physical systems. 

In contrast to the sound wave, the integral over the 
Enskog VACF leads to a finite value: 



knT 



dt{V{t)ViO))^ = 

4-E 



(B13) 



which quantitatively explains the friction we measured 
in our simulations, as discussed in section In con- 
trast to the VACF from sound, the Enskog result cannot 
simply be added to the VACF from Eq. I|B5|I because 
this would violate the t = limit of the VACF, set by 
equipartition. However, because the Enskog contribution 
occurs on times t tcs, where sound and other collective 
modes of (compressible) hydrodynamics are not yet fully 
developed, it may be that Eqs. (jBSp and ljB9|) . derived 
from continuum theory, should be modified on these short 
time-scales. Since the dominant contribution to the in- 
tegral in Eq. ljB7|l is for longer times, the exact form of 
the hydrodynamic contribution IjBSp on these short times 
(t ^ tcs) does not have much influence on the hydrody- 
namic part of the friction. Thus the approximation of 
the short-time VACF by the Enskog result of Eq. IjBlip , 
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and the longer-time {t:$> te, but still t < t^)) VACF by 
Eq. HB5I) may be a fairly reasonable empirical approach. 
This approximation, together with the Green-Kubo rela- 
tion (|B4p . then justifies, a-posteori, the parallel addition 
of frictions in Eq. (|^ . 

Just as was found for the long-time behavior, the short- 
time behavior of the VACF displays important devia- 
tions from the simple Langevin picture. We find that the 
VACF decays much faster, on the order of the sonic time 
rather than the Brownian time Tg. From the tables 
fVl and IVTI one can see how large these differences can be. 
For example, for colloids with a — Ifim, the short-time 
Enskog or sonic decay times are of order 10~^ s, while tb 
is of order 10~^ s, so that the decay of the VACF will be 
dominated by < V{t)V{0) , with its non-exponential 
behavior, over almost the entire time regime. 

An example of the effect of the Enskog contribution on 
the friction or the mean-square displacement is given in 
Fig. CI The integral of the VACF, D{t) from Eq. (EHJ, 
which is directly related to the mean-square displace- 
ment '6 15, is extracted from the simulations described 
in Fig ll3l and compared to 

D"{t)=J^ dt' (^{V{t')V{0)f + {V{t')V{Q)y") (B14) 

calculated using Eqs. (|B5|I and (|B9|I . Both show ex- 
actly the same {t/T^)i scaling for t ^ t^, but the simu- 
lated D{t) is larger because the initial Enskog contribu- 
tion means that D{t) « D^{t) + De, for t ^ te, with 
De = ksT/^E defined in Eq. JBISJ. 

The slow algebraic convergence to the final diffusion 
coefficient, or related friction, also helps explain the scal- 
ing of the finite size effects of Eq. H43|l . A finite box of 
length L means that the measured D{t) is cut off roughly 
at — {L'^ /a?)T^, and because J dtD{t) ~ {T^/t)^ for 
t }t Tv the effect on the diffusion constant can be writ- 
ten as Z) « Do{l — c{L/a)) which explains the scaling 
form of Eq. H43|l . Fig. E]also demonstrates how much 
more rapidly the D{t) from the Langevin Eq. HB3() con- 
verges to its final result D = ksT/^ when compared to 
the simulations and hydrodynamic theories. 

In summary, even though the simple Langevin equa- 
tion (|Bip . without memory, may have pedagogical merit, 
it does a poor job of describing the VACF for colloids in 
solution. 



APPENDIX C: COMPARISON OF SRD TO 
PHYSICAL PARAMETERS 
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FIG. 14: (Color online) Integral D{t) of the normalized VACF 
from Fie. 1181 fsolid line), compared to the integral of Eq. l|B5fl 
(dashed line) , and including the contribution from sound as in 
Eq. I|B14^ (dotted line). The simple Langevin form from the 
integral of Eq. lB3t (dash-dotted line) , is calculated with ^ = 
Q-K-qcjcf in order to have the same asymptotic value as (t) . 
Because in the plot we integrate the normalized VACF, the 
t 00 limit of D{t) for the two theoretical approaches is |io 
in our units, independent of colloid size. The measured D{t) 
has a larger asymptotic value due to the initial effect of the 
Enskog friction. 



Table HVl compares the properties of a pure SRD fluid, 
in the limit of small mean-free path A, to the properties 
of water. 

Table IVl summarizes the main scaling and values for a 
number of different properties of two different size neu- 
trally buoyant colloids in water. To approximate the En- 
skog friction we used a simplified form valid for stick 
boundary conditions in an ideal fluid, as in [58| . which 
ignores the compressibility and other expected complex- 
ities of real water. These predictions should therefore be 
taken with a grain of salt. For simplicity, we have ignored 
any potential Enskog effects on the diffusion constant or 
time-scales. 

Table I VII shows a very similar comparison for two dif- 
ferent size colloids in an SRD fluid. For simplicity we as- 
sume that there are no finite-size effects, which explains 
small differences when compared to similar parameters 
used in the text of the main article. However, in contrast 
to Table we do include the Enskog contribution when 
calculating the total friction and related properties. 



At the end of the day one would like to compare a 
coarse-grained simulation to a real physical colloidal dis- 
persion. To help facilitate such comparisons, we have 
summarized a number of physical parameters and dimen- 
sionless numbers for colloids in water, and also colloids 
in an SRD fluid. 
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TABLE IV: Summary of the physical properties of a 3D SRD fluid in the regime of small reduced mean-free path A. In the 
first column we show how these parameters scale for SRD. To highlight the scaling with the reduced mean-free pathA and the 
number of particles per cell 7, we only show the coUisional contribution to the viscosity in Eqs. @ - 1151 . which dominates for 
small A, and also ignore the factor ( 7, a). In the second column we list the parameter values for our choice of simulation 
parameters in this paper and in refs [52. l58| . and in the last column we compare to the parameter values for II2O at standard 
temperature and pressure. For the SRD the units are in terms of ruf, ao and fesT, so that the time unit is to = ao\/mf /ksT. 
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TABLE V: Physical parameters, dimensionless hydrodynamic 
numbers and time-scales for spherical colloids of radius a — 
0.01/^m (lOnm) and a = 1/im in H2O at standard temperature 
and pressure, moving at a velocity vs ~ lO/xm/s. The mass 
density pc of the colloid is taken to be the same as that of wa- 
ter (the colloid is neutrally buoyant), and the hydrodynamic 
radius a is taken to be the same as the physical radius. Where, 
for notational clarity, the units are not explicitly shown, then 
length is measured in ^m, time in s, and mass in g. 
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TABLE VL Physical parameters, dimensionless hydrody- 
namic numbers and time-scales for a stick boundary colloid 
in an SRD fluid. In the first column, we take the same small 
A limit as in Table HVl to highlight the dominant scaling. In 
the other columns more accurate values for f/vo are used. 
The hydrodynamic numbers and times were estimated with a 
length-scale a ~ a^f ~ o"cc/2. The colloids are neutrally buoy- 
ant, pc = Pf = ymf/ao, and their velocity is vs ~ 0.05ao/to, 
i.e. Ma=0.039. Units as in Table HVI 
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